Análisis exploratorio datos
Isabel Montero Blázquez, Daniel Serrano Serrano, Alejandro Morales Alaminos, Constanza Paz Torres Rivera, Elena Simarro Álvarez y Laura Ruiz Almagro
Primero, cargamos las librerías.
library(ggplot2)
library(ggpubr)
library(summarytools)
library(psych)
library(dplyr)
library(purrr)
library(readxl)
library(scatterplot3d)
library(reshape2)
library(MVN)
library(caret)
library(MASS)
library(biotools)
library(klaR)
library(cluster)
library(ggdendro)
library(gridExtra)
library(tidyverse)
library(factoextra)
library(gridExtra)
library(corrplot)Importamos los datos. Como nuestro excel tiene tres pestañas, lo separamos en 3 data frames distintos. Datos1 para los datos del estudio 1, Datos2 para los del estudio 2 y Datos3 para los del estudio 3.
Datos1<- read_excel("C:\\Users\\Laura\\Desktop\\Datos arreglados.xlsx",1)
#View(Datos1)
Datos2<- read_excel("C:\\Users\\Laura\\Desktop\\Datos arreglados.xlsx",2)
#View(Datos2)
Datos3<- read_excel("C:\\Users\\Laura\\Desktop\\Datos arreglados.xlsx",3)
#View(Datos3)Análisis exploratorio univariante
Datos ausentes (NA). Primero vemos qué datos son los que faltan. Observando las bases de datos, nos damos cuenta de que en el segundo estudio faltan BMC, BMD y longitud de las ratas con ID 222 y 223, y Tb.Sp, Tb.N, BMC,BMD, longitud, GPSurface, GPVolume de las ratas con ID 295 y 298. Tenemos que decidir si quitamos dichas ratas o simplemente omitimos los datos ausentes. Por ahora voy a trabajar omitiendo los datos ausentes. Vamos a realizar un análisis descriptivo separado por estudios.
ESTUDIO 1
Peso inicial
descr(Datos1[,14])## Descriptive Statistics
## Datos1$PesoInicial
## N: 56
##
## PesoInicial
## ----------------- -------------
## Mean 120.13
## Std.Dev 8.75
## Min 102.30
## Q1 113.25
## Median 118.95
## Q3 127.35
## Max 140.80
## MAD 10.30
## IQR 14.05
## CV 0.07
## Skewness 0.22
## SE.Skewness 0.32
## Kurtosis -0.86
## N.Valid 56.00
## Pct.Valid 100.00pesosiniciales1<-Datos1[[14]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosiniciales1, col="pink",xlab="Pesos iniciales",
ylab="Frecuencia", main="Distribución de los pesos iniciales estudio 1")
plot(density(pesosiniciales1), main="Densidad pesos iniciales estudio 1",
xlab="Pesos iniciales", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(pesosiniciales1, main="Pesos iniciales estudio 1", outline=TRUE)
Peso intermedio
descr(Datos1[,15])## Descriptive Statistics
## Datos1$PesoIntermedio
## N: 56
##
## PesoIntermedio
## ----------------- ----------------
## Mean 184.64
## Std.Dev 49.45
## Min 139.20
## Q1 158.65
## Median 167.50
## Q3 181.00
## Max 334.00
## MAD 16.53
## IQR 22.08
## CV 0.27
## Skewness 1.83
## SE.Skewness 0.32
## Kurtosis 2.01
## N.Valid 56.00
## Pct.Valid 100.00pesosintermedios1<-Datos1[[15]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosintermedios1, col="pink", xlab="Pesos intermedios",
ylab="Frecuencia", main="Distribución de los pesos intermedios estudio 1")
plot(density(pesosintermedios1), main="Densidad pesos intermedios estudio 1",
xlab="Pesos intermedios", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(pesosintermedios1, main="Pesos intermedios estudio 1", outline=TRUE)
Peso medición
descr(Datos1[,6])## Descriptive Statistics
## Datos1$PesoFinal
## N: 56
##
## PesoFinal
## ----------------- -----------
## Mean 301.51
## Std.Dev 35.85
## Min 253.20
## Q1 275.80
## Median 294.75
## Q3 314.20
## Max 432.80
## MAD 28.10
## IQR 37.25
## CV 0.12
## Skewness 1.36
## SE.Skewness 0.32
## Kurtosis 2.12
## N.Valid 56.00
## Pct.Valid 100.00pesosmedicion1<-Datos1[[6]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosmedicion1, col="pink", xlab="Pesos medición",
ylab="Frecuencia", main="Distribución de los pesos medición estudio 1")
plot(density(pesosmedicion1), main="Densidad pesos medición estudio 1",
xlab="Pesos medición", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(pesosmedicion1, main="Pesos medición estudio 1", outline=TRUE)
Tb. N
descr(Datos1[,8])## Descriptive Statistics
## Datos1$Tb.N
## N: 56
##
## Tb.N
## ----------------- --------
## Mean 3.14
## Std.Dev 0.57
## Min 1.86
## Q1 2.76
## Median 3.12
## Q3 3.50
## Max 4.62
## MAD 0.54
## IQR 0.73
## CV 0.18
## Skewness 0.08
## SE.Skewness 0.32
## Kurtosis -0.27
## N.Valid 56.00
## Pct.Valid 100.00TbN1<-na.omit(Datos1[[8]])
par(mfrow = c(1, 2),cex=0.7)
hist(TbN1, col="pink", xlab="Número trabecular",
ylab="Frecuencia", main="Distribución del número trabecular estudio 1")
plot(density(TbN1), main="Densidad número trabecular estudio 1",
xlab="Número trabecular", ylab="Densidad") 
par(mfrow = c(1, 1))
boxplot(TbN1, main="Número trabecular estudio 1", outline=TRUE)
BMC
descr(Datos1[,9])## Descriptive Statistics
## Datos1$BMC
## N: 56
##
## BMC
## ----------------- --------
## Mean 0.27
## Std.Dev 0.03
## Min 0.23
## Q1 0.25
## Median 0.27
## Q3 0.29
## Max 0.37
## MAD 0.03
## IQR 0.04
## CV 0.12
## Skewness 1.11
## SE.Skewness 0.32
## Kurtosis 0.67
## N.Valid 56.00
## Pct.Valid 100.00BMC1<-na.omit(Datos1[[9]])
par(mfrow = c(1, 2),cex=0.7)
hist(BMC1, col="pink", xlab="BMC",
ylab="Frecuencia", main="Distribución de BMC estudio 1")
plot(density(BMC1), main="Densidad BMC estudio 1",
xlab="BMC", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(BMC1, main="BMC estudio 1", outline=TRUE)
BMD
descr(Datos1[,10])## Descriptive Statistics
## Datos1$BMD
## N: 56
##
## BMD
## ----------------- --------
## Mean 143.85
## Std.Dev 7.46
## Min 131.95
## Q1 138.50
## Median 142.71
## Q3 147.78
## Max 164.36
## MAD 7.04
## IQR 9.22
## CV 0.05
## Skewness 0.75
## SE.Skewness 0.32
## Kurtosis 0.34
## N.Valid 56.00
## Pct.Valid 100.00BMD1<-na.omit(Datos1[[10]])
par(mfrow = c(1, 2),cex=0.7)
hist(BMD1, col="pink", xlab="BMD",
ylab="Frecuencia", main="Distribución de BMD estudio 1")
plot(density(BMD1), main="BMD estudio 1",
xlab="BMD", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(BMD1, main="BMD estudio 1", outline=TRUE)
Longitud
descr(Datos1[,11])## Descriptive Statistics
## Datos1$Longitud
## N: 56
##
## Longitud
## ----------------- ----------
## Mean 38.92
## Std.Dev 1.10
## Min 37.12
## Q1 38.17
## Median 38.74
## Q3 39.51
## Max 42.09
## MAD 0.93
## IQR 1.31
## CV 0.03
## Skewness 0.69
## SE.Skewness 0.32
## Kurtosis 0.00
## N.Valid 56.00
## Pct.Valid 100.00longitud1<-na.omit(Datos1[[11]])
par(mfrow = c(1, 2),cex=0.7)
hist(longitud1, col="pink", xlab="Longitud",
ylab="Frecuencia", main="Distribución de longitudes estudio 1")
plot(density(longitud1), main="Densidad de la longitud estudio 1",
xlab="Longitud", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(longitud1, main="Longitudes estudio 1", outline=TRUE)
GPSurface
descr(Datos1[,12])## Descriptive Statistics
## Datos1$GPSurface
## N: 56
##
## GPSurface
## ----------------- -----------
## Mean 40.57
## Std.Dev 3.07
## Min 34.36
## Q1 38.78
## Median 40.34
## Q3 42.62
## Max 49.30
## MAD 2.79
## IQR 3.81
## CV 0.08
## Skewness 0.56
## SE.Skewness 0.32
## Kurtosis 0.35
## N.Valid 56.00
## Pct.Valid 100.00GPSurface1<-na.omit(Datos1[[12]])
par(mfrow = c(1, 2),cex=0.7)
hist(GPSurface1, col="pink", xlab="GPSurface",
ylab="Frecuencia", main="Distribución de GPSurface estudio 1")
plot(density(GPSurface1), main="Densidad GPSurface estudio 1",
xlab="GPSurface", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(GPSurface1, main="GPSurface estudio 1", outline=TRUE)
GPVolume
descr(Datos1[,13])## Descriptive Statistics
## Datos1$GPVolume
## N: 56
##
## GPVolume
## ----------------- ----------
## Mean 5.60
## Std.Dev 1.41
## Min 3.35
## Q1 4.66
## Median 5.61
## Q3 6.06
## Max 10.32
## MAD 1.17
## IQR 1.38
## CV 0.25
## Skewness 1.06
## SE.Skewness 0.32
## Kurtosis 1.23
## N.Valid 56.00
## Pct.Valid 100.00GPVolume1<-na.omit(Datos1[[13]])
par(mfrow = c(1, 2),cex=0.7)
hist(GPVolume1, col="pink", xlab="GPVolume",
ylab="Frecuencia", main="Distribución de GPVolume estudio 1")
plot(density(GPVolume1), main="Densidad GPVolume estudio 1",
xlab="GPVolume", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(GPVolume1, main="GPVolume estudio 1", outline=TRUE)
Tb.Sp
descr(Datos1[,7])## Descriptive Statistics
## Datos1$Tb.Sp
## N: 56
##
## Tb.Sp
## ----------------- --------
## Mean 0.33
## Std.Dev 0.07
## Min 0.21
## Q1 0.28
## Median 0.32
## Q3 0.35
## Max 0.55
## MAD 0.05
## IQR 0.07
## CV 0.21
## Skewness 0.93
## SE.Skewness 0.32
## Kurtosis 0.87
## N.Valid 56.00
## Pct.Valid 100.00TbSp1<-na.omit(Datos1[[7]])
par(mfrow = c(1, 2),cex=0.7)
hist(TbSp1, col="pink", xlab="TbSp",
ylab="Frecuencia", main="Distribución del espacio trabecular estudio 1")
plot(density(TbSp1), main="Densidad TbSp estudio 1",
xlab="TbSp", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(TbSp1, main="TbSp estudio 1", outline=TRUE)
Outliers
boxplot(Datos1[,6:15],main="Outliers",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Estandarizamos.
sca<-scale(Datos1[,6:15])
boxplot(sca,main="Outliers",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
ESTUDIO 2
Peso inicial
descr(Datos2[,14])## Descriptive Statistics
## Datos2$PesoInicial
## N: 42
##
## PesoInicial
## ----------------- -------------
## Mean 115.55
## Std.Dev 8.18
## Min 99.40
## Q1 109.40
## Median 116.55
## Q3 119.90
## Max 131.60
## MAD 9.12
## IQR 10.22
## CV 0.07
## Skewness 0.03
## SE.Skewness 0.37
## Kurtosis -0.82
## N.Valid 42.00
## Pct.Valid 100.00pesosiniciales2<-Datos2[[14]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosiniciales2, col="pink", xlab="Pesos iniciales",
ylab="Frecuencia", main="Distribución de los pesos iniciales estudio 2")
plot(density(pesosiniciales2), main="Densidad pesos iniciales estudio 2",
xlab="Pesos iniciales", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(pesosiniciales2, main="Pesos iniciales estudio 2", outline=TRUE)
Peso intermedio
descr(Datos2[,15])## Descriptive Statistics
## Datos2$PesoIntermedio
## N: 42
##
## PesoIntermedio
## ----------------- ----------------
## Mean 201.07
## Std.Dev 92.12
## Min 123.20
## Q1 144.90
## Median 165.60
## Q3 185.50
## Max 397.80
## MAD 29.73
## IQR 39.35
## CV 0.46
## Skewness 1.38
## SE.Skewness 0.37
## Kurtosis 0.15
## N.Valid 42.00
## Pct.Valid 100.00pesosintermedios2<-Datos2[[15]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosintermedios2, col="pink", xlab="Pesos intermedios",
ylab="Frecuencia", main="Distribución de los pesos intermedios estudio 2")
plot(density(pesosintermedios2), main="Densidad pesos intermedios estudio 2",
xlab="Pesos intermedios", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(pesosintermedios2, main="Pesos intermedios estudio 2", outline=TRUE)
Peso medición
descr(Datos2[,6])## Descriptive Statistics
## Datos2$PesoFinal
## N: 42
##
## PesoFinal
## ----------------- -----------
## Mean 358.72
## Std.Dev 54.66
## Min 272.00
## Q1 327.20
## Median 345.30
## Q3 370.60
## Max 480.60
## MAD 36.47
## IQR 42.65
## CV 0.15
## Skewness 0.89
## SE.Skewness 0.37
## Kurtosis -0.24
## N.Valid 42.00
## Pct.Valid 100.00pesosmedicion2<-Datos2[[6]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosmedicion2, col="pink", xlab="Pesos medición",
ylab="Frecuencia", main="Distribución de los pesos medición estudio 2")
plot(density(pesosmedicion2), main="Densidad pesos medición estudio 2",
xlab="Pesos medición", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(pesosmedicion2, main="Pesos medición estudio 2", outline=TRUE)
Tb. N
descr(Datos2[,8])## Descriptive Statistics
## Datos2$Tb.N
## N: 42
##
## Tb.N
## ----------------- -------
## Mean 2.97
## Std.Dev 0.61
## Min 1.67
## Q1 2.58
## Median 3.02
## Q3 3.34
## Max 4.39
## MAD 0.56
## IQR 0.71
## CV 0.20
## Skewness -0.09
## SE.Skewness 0.37
## Kurtosis -0.28
## N.Valid 40.00
## Pct.Valid 95.24TbN2<-na.omit(Datos2[[8]])
par(mfrow = c(1, 2),cex=0.7)
hist(TbN2, col="pink", xlab="Número trabecular",
ylab="Frecuencia", main="Distribución del número trabecular estudio 2")
plot(density(TbN2), main="Densidad número trabecular estudio 2",
xlab="Número trabecular", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(TbN2, main="Número trabecular estudio 2", outline=TRUE)
BMC
descr(Datos2[,9])## Descriptive Statistics
## Datos2$BMC
## N: 42
##
## BMC
## ----------------- -------
## Mean 0.35
## Std.Dev 0.05
## Min 0.26
## Q1 0.32
## Median 0.33
## Q3 0.35
## Max 0.50
## MAD 0.03
## IQR 0.03
## CV 0.15
## Skewness 1.19
## SE.Skewness 0.38
## Kurtosis 0.72
## N.Valid 38.00
## Pct.Valid 90.48BMC2<-na.omit(Datos2[[9]])
par(mfrow = c(1, 2),cex=0.7)
hist(BMC2, col="pink", xlab="BMC",
ylab="Frecuencia", main="Distribución de BMC estudio 2")
plot(density(BMC2), main="Densidad BMC estudio 2",
xlab="BMC", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(BMC2, main="BMC estudio 2", outline=TRUE)
BMD
descr(Datos2[,10])## Descriptive Statistics
## Datos2$BMD
## N: 42
##
## BMD
## ----------------- --------
## Mean 159.34
## Std.Dev 12.07
## Min 137.87
## Q1 151.25
## Median 157.38
## Q3 166.30
## Max 186.94
## MAD 9.45
## IQR 14.53
## CV 0.08
## Skewness 0.66
## SE.Skewness 0.38
## Kurtosis -0.36
## N.Valid 38.00
## Pct.Valid 90.48BMD2<-na.omit(Datos2[[10]])
par(mfrow = c(1, 2),cex=0.7)
hist(BMD2, col="pink", xlab="BMD",
ylab="Frecuencia", main="Distribución de BMD estudio 2")
plot(density(BMD2), main="BMD estudio 2",
xlab="BMD", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(BMD2, main="BMD estudio 2", outline=TRUE)
Longitud
descr(Datos2[,11])## Descriptive Statistics
## Datos2$Longitud
## N: 42
##
## Longitud
## ----------------- ----------
## Mean 42.04
## Std.Dev 1.47
## Min 39.29
## Q1 41.32
## Median 41.74
## Q3 42.73
## Max 45.52
## MAD 0.73
## IQR 1.27
## CV 0.04
## Skewness 0.65
## SE.Skewness 0.38
## Kurtosis 0.04
## N.Valid 38.00
## Pct.Valid 90.48longitud2<-na.omit(Datos2[[11]])
par(mfrow = c(1, 2),cex=0.7)
hist(longitud2, col="pink", xlab="Longitud",
ylab="Frecuencia", main="Distribución de longitudes estudio 2")
plot(density(longitud2), main="Densidad de la longitud estudio 2",
xlab="Longitud", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(longitud2, main="Longitudes estudio 2", outline=TRUE)
GPSurface
descr(Datos2[,12])## Descriptive Statistics
## Datos2$GPSurface
## N: 42
##
## GPSurface
## ----------------- -----------
## Mean 44.12
## Std.Dev 4.44
## Min 35.71
## Q1 41.27
## Median 44.03
## Q3 45.63
## Max 54.68
## MAD 3.62
## IQR 4.28
## CV 0.10
## Skewness 0.46
## SE.Skewness 0.37
## Kurtosis -0.15
## N.Valid 40.00
## Pct.Valid 95.24GPSurface2<-na.omit(Datos2[[12]])
par(mfrow = c(1, 2),cex=0.7)
hist(GPSurface2, col="pink", xlab="GPSurface",
ylab="Frecuencia", main="Distribución de GPSurface estudio 2")
plot(density(GPSurface2), main="Densidad GPSurface estudio 2",
xlab="GPSurface", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(GPSurface2, main="GPSurface estudio 2", outline=TRUE)
GPVolume
descr(Datos2[,13])## Descriptive Statistics
## Datos2$GPVolume
## N: 42
##
## GPVolume
## ----------------- ----------
## Mean 5.57
## Std.Dev 0.97
## Min 3.59
## Q1 4.79
## Median 5.52
## Q3 6.15
## Max 7.85
## MAD 0.99
## IQR 1.33
## CV 0.17
## Skewness 0.26
## SE.Skewness 0.37
## Kurtosis -0.41
## N.Valid 40.00
## Pct.Valid 95.24GPVolume2<-na.omit(Datos2[[13]])
par(mfrow = c(1, 2),cex=0.7)
hist(GPVolume2, col="pink", xlab="GPVolume",
ylab="Frecuencia", main="Distribución de GPVolume estudio 2")
plot(density(GPVolume2), main="Densidad GPVolume estudio 2",
xlab="GPVolume", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(GPVolume2, main="GPVolume estudio 2", outline=TRUE)
Tb.Sp
descr(Datos2[,7])## Descriptive Statistics
## Datos2$Tb.Sp
## N: 42
##
## Tb.Sp
## ----------------- -------
## Mean 0.35
## Std.Dev 0.09
## Min 0.21
## Q1 0.29
## Median 0.33
## Q3 0.38
## Max 0.60
## MAD 0.06
## IQR 0.08
## CV 0.25
## Skewness 1.09
## SE.Skewness 0.37
## Kurtosis 0.80
## N.Valid 40.00
## Pct.Valid 95.24TbSp2<-na.omit(Datos2[[7]])
par(mfrow = c(1, 2),cex=0.7)
hist(TbSp2, col="pink", xlab="TbSp",
ylab="Frecuencia", main="Distribución del espacio trabecular estudio 2")
plot(density(TbSp2), main="Densidad TbSp estudio 2",
xlab="TbSp", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(TbSp2, main="TbSp estudio 2", outline=TRUE)
Outliers
boxplot(Datos2[,6:15],main="Outliers",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Estandarizamos.
sca<-scale(Datos2[,6:15])
boxplot(sca,main="Outliers",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
ESTUDIO 3
Peso inicial
descr(Datos3[,14])## Descriptive Statistics
## Datos3$PesoInicial
## N: 36
##
## PesoInicial
## ----------------- -------------
## Mean 71.76
## Std.Dev 4.56
## Min 63.00
## Q1 69.35
## Median 71.65
## Q3 73.30
## Max 86.80
## MAD 2.74
## IQR 3.83
## CV 0.06
## Skewness 0.95
## SE.Skewness 0.39
## Kurtosis 2.16
## N.Valid 36.00
## Pct.Valid 100.00pesosiniciales3<-Datos3[[14]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosiniciales3, col="pink", xlab="Pesos iniciales",
ylab="Frecuencia", main="Distribución de los pesos iniciales estudio 3")
plot(density(pesosiniciales3), main="Densidad pesos iniciales estudio 3",
xlab="Pesos iniciales", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(pesosiniciales3, main="Pesos iniciales estudio 3", outline=TRUE)
Peso intermedio
descr(Datos3[,15])## Descriptive Statistics
## Datos3$PesoIntermedio
## N: 36
##
## PesoIntermedio
## ----------------- ----------------
## Mean 127.51
## Std.Dev 51.49
## Min 77.00
## Q1 95.05
## Median 104.25
## Q3 166.65
## Max 234.00
## MAD 17.27
## IQR 54.10
## CV 0.40
## Skewness 1.04
## SE.Skewness 0.39
## Kurtosis -0.73
## N.Valid 36.00
## Pct.Valid 100.00pesosintermedios3<-Datos3[[15]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosintermedios3, col="pink", xlab="Pesos intermedios",
ylab="Frecuencia", main="Distribución de los pesos intermedios estudio 3")
plot(density(pesosintermedios3), main="Densidad pesos intermedios estudio 3",
xlab="Pesos intermedios", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(pesosintermedios3, main="Pesos intermedios estudio 3", outline=TRUE)
Peso medición
descr(Datos3[,6])## Descriptive Statistics
## Datos3$PesoFinal
## N: 36
##
## PesoFinal
## ----------------- -----------
## Mean 238.28
## Std.Dev 47.10
## Min 163.00
## Q1 208.00
## Median 223.00
## Q3 275.65
## Max 330.00
## MAD 26.02
## IQR 61.32
## CV 0.20
## Skewness 0.68
## SE.Skewness 0.39
## Kurtosis -0.84
## N.Valid 36.00
## Pct.Valid 100.00pesosmedicion3<-Datos3[[6]]
par(mfrow = c(1, 2),cex=0.7)
hist(pesosmedicion3, col="pink", xlab="Pesos medición",
ylab="Frecuencia", main="Distribución de los pesos medición estudio 3")
plot(density(pesosmedicion3), main="Densidad pesos medición estudio 3",
xlab="Pesos medición", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(pesosmedicion3, main="Pesos medición estudio 3", outline=TRUE)
Tb. N
descr(Datos3[,8])## Descriptive Statistics
## Datos3$Tb.N
## N: 36
##
## Tb.N
## ----------------- --------
## Mean 3.29
## Std.Dev 0.61
## Min 1.79
## Q1 2.78
## Median 3.48
## Q3 3.75
## Max 4.24
## MAD 0.69
## IQR 0.95
## CV 0.18
## Skewness -0.45
## SE.Skewness 0.39
## Kurtosis -0.73
## N.Valid 36.00
## Pct.Valid 100.00TbN3<-na.omit(Datos3[[8]])
par(mfrow = c(1, 2),cex=0.7)
hist(TbN3, col="pink", xlab="Número trabecular",
ylab="Frecuencia", main="Distribución del número trabecular estudio 3")
plot(density(TbN3), main="Densidad número trabecular estudio 3",
xlab="Número trabecular", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(TbN3, main="Número trabecular estudio 3", outline=TRUE)
BMC
descr(Datos3[,9])## Descriptive Statistics
## Datos3$BMC
## N: 36
##
## BMC
## ----------------- --------
## Mean 0.20
## Std.Dev 0.04
## Min 0.15
## Q1 0.17
## Median 0.19
## Q3 0.24
## Max 0.30
## MAD 0.03
## IQR 0.06
## CV 0.21
## Skewness 0.76
## SE.Skewness 0.39
## Kurtosis -0.73
## N.Valid 36.00
## Pct.Valid 100.00BMC3<-na.omit(Datos3[[9]])
par(mfrow = c(1, 2),cex=0.7)
hist(BMC3, col="pink", xlab="BMC",
ylab="Frecuencia", main="Distribución de BMC estudio 3")
plot(density(BMC3), main="Densidad BMC estudio 3",
xlab="BMC", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(BMC3, main="BMC estudio 3", outline=TRUE)
BMD
descr(Datos3[,10])## Descriptive Statistics
## Datos3$BMD
## N: 36
##
## BMD
## ----------------- --------
## Mean 123.57
## Std.Dev 11.05
## Min 105.62
## Q1 113.55
## Median 122.86
## Q3 131.92
## Max 146.00
## MAD 13.56
## IQR 18.03
## CV 0.09
## Skewness 0.42
## SE.Skewness 0.39
## Kurtosis -0.93
## N.Valid 36.00
## Pct.Valid 100.00BMD3<-na.omit(Datos3[[10]])
par(mfrow = c(1, 2),cex=0.7)
hist(BMD3, col="pink", xlab="BMD",
ylab="Frecuencia", main="Distribución de BMD estudio 3")
plot(density(BMD3), main="BMD estudio 3",
xlab="BMD", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(BMD3, main="BMD estudio 3", outline=TRUE)
Longitud
descr(Datos3[,11])## Descriptive Statistics
## Datos3$Longitud
## N: 36
##
## Longitud
## ----------------- ----------
## Mean 34.84
## Std.Dev 1.59
## Min 32.45
## Q1 33.61
## Median 34.30
## Q3 36.16
## Max 37.89
## MAD 1.19
## IQR 2.46
## CV 0.05
## Skewness 0.55
## SE.Skewness 0.39
## Kurtosis -1.02
## N.Valid 36.00
## Pct.Valid 100.00longitud3<-na.omit(Datos3[[11]])
par(mfrow = c(1, 2),cex=0.7)
hist(longitud3, col="pink", xlab="Longitud",
ylab="Frecuencia", main="Distribución de longitudes estudio 3")
plot(density(longitud3), main="Densidad de la longitud estudio 3",
xlab="Longitud", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(longitud3, main="Longitudes estudio 3", outline=TRUE)
GPSurface
descr(Datos3[,12])## Descriptive Statistics
## Datos3$GPSurface
## N: 36
##
## GPSurface
## ----------------- -----------
## Mean 35.90
## Std.Dev 3.71
## Min 28.96
## Q1 33.73
## Median 34.47
## Q3 38.31
## Max 44.20
## MAD 1.93
## IQR 4.43
## CV 0.10
## Skewness 0.69
## SE.Skewness 0.39
## Kurtosis -0.36
## N.Valid 36.00
## Pct.Valid 100.00GPSurface3<-na.omit(Datos3[[12]])
par(mfrow = c(1, 2),cex=0.7)
hist(GPSurface3, col="pink", xlab="GPSurface",
ylab="Frecuencia", main="Distribución de GPSurface estudio 3")
plot(density(GPSurface3), main="Densidad GPSurface estudio 3",
xlab="GPSurface", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(GPSurface3, main="GPSurface estudio 3", outline=TRUE)
GPVolume
descr(Datos3[,13])## Descriptive Statistics
## Datos3$GPVolume
## N: 36
##
## GPVolume
## ----------------- ----------
## Mean 5.65
## Std.Dev 0.86
## Min 4.04
## Q1 4.99
## Median 5.57
## Q3 6.32
## Max 7.58
## MAD 0.93
## IQR 1.28
## CV 0.15
## Skewness 0.13
## SE.Skewness 0.39
## Kurtosis -0.69
## N.Valid 36.00
## Pct.Valid 100.00GPVolume3<-na.omit(Datos3[[13]])
par(mfrow = c(1, 2),cex=0.7)
hist(GPVolume3, col="pink", xlab="GPVolume",
ylab="Frecuencia", main="Distribución de GPVolume estudio 3")
plot(density(GPVolume3), main="Densidad GPVolume estudio 3",
xlab="GPVolume", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(GPVolume3, main="GPVolume estudio 3", outline=TRUE)
Tb.Sp
descr(Datos3[,7])## Descriptive Statistics
## Datos3$Tb.Sp
## N: 36
##
## Tb.Sp
## ----------------- --------
## Mean 0.31
## Std.Dev 0.07
## Min 0.22
## Q1 0.26
## Median 0.28
## Q3 0.36
## Max 0.57
## MAD 0.05
## IQR 0.10
## CV 0.24
## Skewness 1.25
## SE.Skewness 0.39
## Kurtosis 1.57
## N.Valid 36.00
## Pct.Valid 100.00TbSp3<-na.omit(Datos3[[7]])
par(mfrow = c(1, 2),cex=0.7)
hist(TbSp3, col="pink", xlab="TbSp",
ylab="Frecuencia", main="Distribución del espacio trabecular estudio 3")
plot(density(TbSp3), main="Densidad TbSp estudio 3",
xlab="TbSp", ylab="Densidad")
par(mfrow = c(1, 1))
boxplot(TbSp3, main="TbSp estudio 3", outline=TRUE)
Outliers
boxplot(Datos3[,6:15],main="Outliers",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Estandarizamos.
sca<-scale(Datos3[,6:15])
boxplot(sca,main="Outliers",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Vamos a superponer los histogramas para ver cómo cambian las variables en función de la dosis.
Estudio 1
par(mfrow = c(1, 2),cex=0.7)
p1_1<-ggplot(data=Datos1,aes(x=PesoInicial,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p2_1<-ggplot(data=Datos1,aes(x=PesoIntermedio,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p3_1<-ggplot(data=Datos1,aes(x=PesoFinal,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p4_1<-ggplot(data=Datos1,aes(x=Tb.N,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p5_1<-ggplot(data=Datos1,aes(x=Tb.Sp,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p6_1<-ggplot(data=Datos1,aes(x=BMC,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p7_1<-ggplot(data=Datos1,aes(x=BMD,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p8_1<-ggplot(data=Datos1,aes(x=Longitud,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p9_1<-ggplot(data=Datos1,aes(x=GPSurface,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p10_1<-ggplot(data=Datos1,aes(x=GPVolume,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
ggarrange(p1_1,p2_1, p3_1,p4_1,p5_1,p6_1,p7_1,p8_1,p9_1,p10_1, nrow = 4, ncol=3, common.legend = TRUE)
Estudio 2
Datos2.2<-na.omit(Datos2)
par(mfrow = c(1, 2),cex=0.7)
p1_2<-ggplot(data=Datos2.2,aes(x=PesoInicial,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p2_2<-ggplot(data=Datos2.2,aes(x=PesoIntermedio,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p3_2<-ggplot(data=Datos2.2,aes(x=PesoFinal,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p4_2<-ggplot(data=Datos2.2,aes(x=Tb.N,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p5_2<-ggplot(data=Datos2.2,aes(x=Tb.Sp,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p6_2<-ggplot(data=Datos2.2,aes(x=BMC,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p7_2<-ggplot(data=Datos2.2,aes(x=BMD,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p8_2<-ggplot(data=Datos2.2,aes(x=Longitud,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p9_2<-ggplot(data=Datos2.2,aes(x=GPSurface,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p10_2<-ggplot(data=Datos2.2,aes(x=GPVolume,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
ggarrange(p1_2,p2_2, p3_2,p4_2,p5_2,p6_2,p7_2,p8_2,p9_2,p10_2, nrow = 4, ncol=3, common.legend = TRUE)
Estudio 3
par(mfrow = c(1, 2),cex=0.7)
p1_3<-ggplot(data=Datos3,aes(x=PesoInicial,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p2_3<-ggplot(data=Datos3,aes(x=PesoIntermedio,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p3_3<-ggplot(data=Datos3,aes(x=PesoFinal,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p4_3<-ggplot(data=Datos3,aes(x=Tb.N,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p5_3<-ggplot(data=Datos3,aes(x=Tb.Sp,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p6_3<-ggplot(data=Datos3,aes(x=BMC,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p7_3<-ggplot(data=Datos3,aes(x=BMD,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p8_3<-ggplot(data=Datos3,aes(x=Longitud,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p9_3<-ggplot(data=Datos3,aes(x=GPSurface,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
p10_3<-ggplot(data=Datos3,aes(x=GPVolume,fill=factor(Dosis))) +
geom_histogram(position = "identity", alpha = 0.5, bins=20)
ggarrange(p1_3,p2_3, p3_3,p4_3,p5_3,p6_3,p7_3,p8_3,p9_3,p10_3, nrow = 4, ncol=3, common.legend = TRUE)
En general parece que el peso intermedio, el BMC y la longitud separan mejor al grupo de control. Ahora vamos a realizar un análisis exploratorio visual bivariante.
Análisis exploratorio bivariante
Estudio 1
pairs(Datos1[, c("PesoInicial", "PesoIntermedio", "PesoFinal", "Tb.N", "Tb.Sp", "BMC", "BMD", "Longitud", "GPSurface", "GPVolume")],
col = rainbow(length(levels(factor(Datos1$Dosis))))[factor(Datos1$Dosis)],
pch = 19)
legend("topright",
legend = levels(factor(Datos1$Dosis)),
fill = rainbow(length(levels(factor(Datos1$Dosis)))))
Correlaciones
par(mar = c(3, 1, 4, 1))
corrplot(cor(Datos1[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Estudio 1", side = 3, line = 3, adj = 0.3, cex = 1.1)
Correlaciones por dosis
dato1dosis0<-subset(Datos1,Dosis==0)
dato1dosis1<-subset(Datos1,Dosis==1)
dato1dosis8<-subset(Datos1,Dosis==8)
dato1dosis10<-subset(Datos1,Dosis==10)
dato1dosis80<-subset(Datos1,Dosis==80)
dato1dosis100<-subset(Datos1,Dosis==100)
par(mar = c(3, 1, 4, 1))
corrplot(cor(dato1dosis0[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 0", side = 3, line = 3, adj = 0.3, cex = 1.1)
corrplot(cor(dato1dosis1[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 1", side = 3, line = 3, adj = 0.3, cex = 1.1)
corrplot(cor(dato1dosis8[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 8", side = 3, line = 3, adj = 0.3, cex = 1.1)
corrplot(cor(dato1dosis10[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 10", side = 3, line = 3, adj = 0.3, cex = 1.1)
corrplot(cor(dato1dosis80[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 80", side = 3, line = 3, adj = 0.3, cex = 1.1)
corrplot(cor(dato1dosis100[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 100", side = 3, line = 3, adj = 0.3, cex = 1.1)
Estudio 2
Omitimos los valores perdidos del Estudio 2 para hacer la gráfica. No sería un error eliminarlos ya que son menos del 5%.
Datos2.2<-na.omit(Datos2)
pairs(Datos2.2[, c("PesoInicial", "PesoIntermedio", "PesoFinal", "Tb.N", "Tb.Sp", "BMC", "BMD", "Longitud", "GPSurface", "GPVolume")],
col = rainbow(length(levels(factor(Datos2.2$Dosis))))[factor(Datos2.2$Dosis)],pch = 19)
legend("topright",
legend = levels(factor(Datos2.2$Dosis)),
fill = rainbow(length(levels(factor(Datos2.2$Dosis)))))
Correlaciones
par(mar = c(3, 1, 4, 1))
corrplot(cor(Datos2.2[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Estudio 2", side = 3, line = 3, adj = 0.3, cex = 1.1)
Correlaciones por dosis
dato2dosis0<-subset(Datos2.2,Dosis==0)
dato2dosis80<-subset(Datos2.2,Dosis==80)
dato2dosis400<-subset(Datos2.2,Dosis==400)
par(mar = c(3, 1, 4, 1))
corrplot(cor(dato2dosis0[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 0", side = 3, line = 3, adj = 0.3, cex = 1.1)
corrplot(cor(dato2dosis80[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 80", side = 3, line = 3, adj = 0.3, cex = 1.1)
corrplot(cor(dato2dosis400[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Dosis 400", side = 3, line = 3, adj = 0.3, cex = 1.1)
Estudio 3
pairs(Datos3[, c("PesoInicial", "PesoIntermedio", "PesoFinal", "Tb.N", "Tb.Sp", "BMC", "BMD", "Longitud", "GPSurface", "GPVolume")],
col = rainbow(length(levels(factor(Datos3$Dosis))))[factor(Datos3$Dosis)],pch = 19)
legend("topright",
legend = levels(factor(Datos3$Dosis)),
fill = rainbow(length(levels(factor(Datos3$Dosis)))))
#### Correlaciones
par(mar = c(3, 1, 4, 1))
corrplot(cor(Datos3[6:15]), order = "hclust", tl.col='black', tl.cex=0.6)
mtext("Estudio 3", side = 3, line = 2, adj = 0.3, cex = 1.1)
Correlaciones por dosis
dato3dosis0<-subset(Datos3,Dosis==0)
dato3dosis500<-subset(Datos3,Dosis==500)
dato3dosis1000<-subset(Datos3,Dosis==1000)
par(mar = c(3, 1, 4, 1))
corrplot(cor(dato3dosis0[6:15]), order = "hclust", tl.col='black', tl.cex=0.5)
mtext("Dosis 0", side = 3, line = 3, adj = 0.3, cex = 1.1)
corrplot(cor(dato3dosis500[6:15]), order = "hclust", tl.col='black', tl.cex=0.5)
mtext("Dosis 500", side = 3, line = 3, adj = 0.3, cex = 1.1)
corrplot(cor(dato3dosis1000[6:15]), order = "hclust", tl.col='black', tl.cex=0.5)
mtext("Dosis 1000", side = 3, line = 3, adj = 0.3, cex = 1.1)
Análisis exploratorio trivariante
Finalmente, nos queda hacer un análisis exploratorio visual trivariante. Comparamos peso intermedio, longitud, BMC, BMD, GPSurface y Tb.Sp que en los gráficos de histogramas parecían los que más información nos aportaban.
Estudio 1
Datos1Factor <- Datos1
Datos2.2Factor <- Datos2.2
Datos3Factor <- Datos3
Datos1Factor$Dosis <- factor(Datos1Factor$Dosis)
Datos2.2Factor$Dosis <- factor(Datos2.2Factor$Dosis)
Datos3Factor$Dosis <- factor(Datos3Factor$Dosis)
par(mfrow = c(1, 2),cex=0.7)
#####
scatterplot3d(Datos1$PesoInicial, Datos1$PesoIntermedio, Datos1$PesoFinal,
color = rainbow(length(levels(Datos1Factor$Dosis)))[as.factor(Datos1$Dosis)],
pch = 19, grid = TRUE, xlab = "Peso Inicial", ylab = "Peso Intermedio", zlab = "Peso Final",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(factor(Datos1$Dosis)), fill = rainbow(length(levels(factor(Datos1$Dosis)))))
#####
scatterplot3d(Datos1$Longitud, Datos1$PesoIntermedio, Datos1$BMC,
color = rainbow(length(levels(Datos1Factor$Dosis)))[as.factor(Datos1$Dosis)], pch = 19,
grid = TRUE, xlab = "Longitud", ylab = "Peso Intermedio", zlab = "BMC",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))
####
scatterplot3d(Datos1Factor$GPSurface, Datos1Factor$PesoIntermedio, Datos1Factor$BMC,
color = rainbow(length(levels(Datos1Factor$Dosis)))[Datos1Factor$Dosis], pch = 19,
grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "BMC",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))
####
scatterplot3d(Datos1Factor$GPSurface, Datos1Factor$PesoIntermedio, Datos1Factor$Longitud,
color = rainbow(length(levels(Datos1Factor$Dosis)))[Datos1Factor$Dosis], pch = 19,
grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "Longitud",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))
####
scatterplot3d(Datos1Factor$BMD, Datos1Factor$PesoIntermedio, Datos1Factor$Tb.Sp,
color = rainbow(length(levels(Datos1Factor$Dosis)))[Datos1Factor$Dosis], pch = 19,
grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "Tb.Sp",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))
#####
scatterplot3d(Datos1Factor$BMD, Datos1Factor$PesoIntermedio, Datos1Factor$Longitud,
color = rainbow(length(levels(Datos1Factor$Dosis)))[Datos1Factor$Dosis], pch = 19,
grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "Longitud",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))
#####
scatterplot3d(Datos1Factor$BMC, Datos1Factor$PesoIntermedio, Datos1Factor$Tb.Sp,
color = rainbow(length(levels(Datos1Factor$Dosis)))[Datos1Factor$Dosis], pch = 19,
grid = TRUE, xlab = "BMC", ylab = "Peso Intermedio", zlab = "Tb.Sp",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))
#####
scatterplot3d(Datos1Factor$BMD, Datos1Factor$PesoIntermedio, Datos1Factor$GPSurface,
color = rainbow(length(levels(Datos1Factor$Dosis)))[Datos1Factor$Dosis], pch = 19,
grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "GPSurface",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos1Factor$Dosis), fill = rainbow(length(levels(Datos1Factor$Dosis))))
Estudio 2
par(mfrow = c(1, 2),cex=0.7)
scatterplot3d(Datos2.2$PesoInicial, Datos2.2$PesoIntermedio, Datos2.2$PesoFinal,
color = rainbow(length(levels(factor(Datos2.2$Dosis))))[as.factor(Datos2.2$Dosis)],
pch = 19, grid = TRUE, xlab = "Peso Inicial", ylab = "Peso Intermedio", zlab = "Peso Final",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(factor(Datos2.2$Dosis)), fill = rainbow(length(levels(factor(Datos2.2$Dosis)))))
#####
scatterplot3d(Datos2.2Factor$Longitud, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$BMC,
color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
grid = TRUE, xlab = "Longitud", ylab = "Peso Intermedio", zlab = "BMC",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))
#####
scatterplot3d(Datos2.2Factor$GPSurface, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$BMC,
color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "BMC",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))
######
scatterplot3d(Datos2.2Factor$GPSurface, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$Longitud,
color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "Longitud",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))
#######
scatterplot3d(Datos2.2Factor$BMD, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$Tb.Sp,
color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "Tb.Sp",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))
######
scatterplot3d(Datos2.2Factor$BMD, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$Longitud,
color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "Longitud",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))
#######
scatterplot3d(Datos2.2Factor$BMC, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$Tb.Sp,
color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
grid = TRUE, xlab = "BMC", ylab = "Peso Intermedio", zlab = "Tb.Sp",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))
######
scatterplot3d(Datos2.2Factor$BMD, Datos2.2Factor$PesoIntermedio, Datos2.2Factor$GPSurface,
color = rainbow(length(levels(Datos2.2Factor$Dosis)))[Datos2.2Factor$Dosis], pch = 19,
grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "GPSurface",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos2.2Factor$Dosis), fill = rainbow(length(levels(Datos2.2Factor$Dosis))))
Estudio 3
par(mfrow = c(1, 2),cex=0.7)
scatterplot3d(Datos3$PesoInicial, Datos3$PesoIntermedio, Datos3$PesoFinal,
color = rainbow(length(levels(factor(Datos3$Dosis))))[as.factor(Datos3$Dosis)],
pch = 19, grid = TRUE, xlab = "Peso Inicial", ylab = "Peso Intermedio", zlab = "Peso Final",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(factor(Datos3$Dosis)), fill = rainbow(length(levels(factor(Datos3$Dosis)))))
#####
scatterplot3d(Datos3Factor$Longitud, Datos3Factor$PesoIntermedio, Datos3Factor$BMC,
color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
grid = TRUE, xlab = "Longitud", ylab = "Peso Intermedio", zlab = "BMC",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))
######
scatterplot3d(Datos3Factor$GPSurface, Datos3Factor$PesoIntermedio, Datos3Factor$BMC,
color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "BMC",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))
#######
scatterplot3d(Datos3Factor$GPSurface, Datos3Factor$PesoIntermedio, Datos3Factor$Longitud,
color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "Longitud",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))
#########
scatterplot3d(Datos3Factor$BMD, Datos3Factor$PesoIntermedio, Datos3Factor$Tb.Sp,
color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "Tb.Sp",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))
#########
scatterplot3d(Datos3Factor$BMD, Datos3Factor$PesoIntermedio, Datos3Factor$Longitud,
color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
grid = TRUE, xlab = "BMD", ylab = "Peso Intermedio", zlab = "Longitud",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))
#######
scatterplot3d(Datos3Factor$BMD, Datos3Factor$PesoIntermedio, Datos3Factor$Tb.Sp,
color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
grid = TRUE, xlab = "GPSurface", ylab = "Peso Intermedio", zlab = "Tb.Sp",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))
#######
scatterplot3d(Datos3Factor$BMC, Datos3Factor$PesoIntermedio, Datos3Factor$GPSurface,
color = rainbow(length(levels(Datos3Factor$Dosis)))[Datos3Factor$Dosis], pch = 19,
grid = TRUE, xlab = "BMC", ylab = "Peso Intermedio", zlab = "GPSurface",
angle = 65, cex.axis = 0.6)
legend("topright", bty = "n", cex = 0.8, title = "Dosis",
legend = levels(Datos3Factor$Dosis), fill = rainbow(length(levels(Datos3Factor$Dosis))))
Tratamiento de outliers
Todo el análisis que hemos hecho anteriormente ha sido sin quitar los outliers. Hay tres formas de gestionar los outliers: dejarlos, eliminarlos y cambiarlos por la media.
Estudio 1
Eliminar los outliers
# Función para eliminar outliers de columnas específicas
eliminar_outliers_columnas <- function(datos, columnas) {
# Inicializar el dataframe sin outliers
datos_sin_outliers <- datos
for (columna in columnas) {
# Calcular los cuartiles y el IQR de la columna
Q1 <- quantile(datos[[columna]], 0.25)
Q3 <- quantile(datos[[columna]], 0.75)
IQR <- Q3 - Q1
# Definir los límites inferior y superior
limite_inferior <- Q1 - 1.5 * IQR
limite_superior <- Q3 + 1.5 * IQR
# Filtrar las filas que están dentro de los límites
datos_sin_outliers <- datos_sin_outliers[datos_sin_outliers[[columna]] >= limite_inferior & datos_sin_outliers[[columna]] <= limite_superior, ]
}
return(datos_sin_outliers)
}
columnas_a_usar <- c(6:15)
Datos1_sin_outliers <- eliminar_outliers_columnas(Datos1, columnas_a_usar)
#print(Datos1_sin_outliers)Comprobemos gráficamente la eliminación de los outliers.
par(mfrow = c(1, 2))
boxplot(Datos1[,6:15],main="Outliers antes de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
boxplot(Datos1_sin_outliers[,6:15],main="Outliers después de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Estandarizamos.
par(mfrow = c(1, 2))
sca<-scale(Datos1[,6:15])
boxplot(sca,main="Outliers antes de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
sca_sin_outliers<-scale(Datos1_sin_outliers[,6:15])
boxplot(sca_sin_outliers,main="Outliers después de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Cambiar los outliers por la media
outlier<-function(data,na.rm=T){
H<-1.5*IQR(data)
data[data<quantile(data,0.25,na.rm = T)-H]<-NA
data[data>quantile(data,0.75, na.rm = T)+H]<-NA
data[is.na(data)]<-mean(data, na.rm = T)
H<-1.5*IQR(data)
if (TRUE %in% (data<quantile(data,0.25,na.rm = T)-H) |
TRUE %in% (data>quantile(data,0.75,na.rm = T)+H))
outlier(data)
else
return(data)
}
Datos1_con_media<-Datos1
Datos1_con_media[[6]]<-outlier(Datos1_con_media[[6]])
Datos1_con_media[[7]]<-outlier(Datos1_con_media[[7]])
Datos1_con_media[[8]]<-outlier(Datos1_con_media[[8]])
Datos1_con_media[[9]]<-outlier(Datos1_con_media[[9]])
Datos1_con_media[[10]]<-outlier(Datos1_con_media[[10]])
Datos1_con_media[[11]]<-outlier(Datos1_con_media[[11]])
Datos1_con_media[[12]]<-outlier(Datos1_con_media[[12]])
Datos1_con_media[[13]]<-outlier(Datos1_con_media[[13]])
Datos1_con_media[[14]]<-outlier(Datos1_con_media[[14]])
Datos1_con_media[[15]]<-outlier(Datos1_con_media[[15]]) Comprobemos gráficamente el cambio de los outliers por la media.
par(mfrow = c(1, 2))
boxplot(Datos1[,6:15],main="Outliers antes de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
boxplot(Datos1_con_media[,6:15],main="Outliers después de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Estandarizamos.
par(mfrow = c(1, 2))
sca<-scale(Datos1[,6:15])
boxplot(sca,main="Outliers antes de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
sca_sin_outliers1<-scale(Datos1_con_media[,6:15])
boxplot(sca_sin_outliers,main="Outliers después de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Estudio 2
Eliminar los outliers
# Función para eliminar outliers de columnas específicas
eliminar_outliers_columnas <- function(datos, columnas) {
# Inicializar el dataframe sin outliers
datos_sin_outliers <- datos
for (columna in columnas) {
# Calcular los cuartiles y el IQR de la columna
Q1 <- quantile(datos[[columna]], 0.25)
Q3 <- quantile(datos[[columna]], 0.75)
IQR <- Q3 - Q1
# Definir los límites inferior y superior
limite_inferior <- Q1 - 1.5 * IQR
limite_superior <- Q3 + 1.5 * IQR
# Filtrar las filas que están dentro de los límites
datos_sin_outliers <- datos_sin_outliers[datos_sin_outliers[[columna]] >= limite_inferior & datos_sin_outliers[[columna]] <= limite_superior, ]
}
return(datos_sin_outliers)
}
columnas_a_usar <- c(6:15)
Datos2.2_sin_outliers <- eliminar_outliers_columnas(Datos2.2, columnas_a_usar)
#print(Datos2.2_sin_outliers)Comprobemos gráficamente la eliminación de los outliers.
par(mfrow = c(1, 2))
boxplot(Datos2.2[,6:15],main="Outliers antes de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
boxplot(Datos2.2_sin_outliers[,6:15],main="Outliers después de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Estandarizamos.
par(mfrow = c(1, 2))
sca<-scale(Datos2.2[,6:15])
boxplot(sca,main="Outliers antes de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
sca_sin_outliers<-scale(Datos2.2_sin_outliers[,6:15])
boxplot(sca_sin_outliers,main="Outliers después de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Cambiar los outliers por la media
outlier<-function(data,na.rm=T){
H<-1.5*IQR(data)
data[data<quantile(data,0.25,na.rm = T)-H]<-NA
data[data>quantile(data,0.75, na.rm = T)+H]<-NA
data[is.na(data)]<-mean(data, na.rm = T)
H<-1.5*IQR(data)
if (TRUE %in% (data<quantile(data,0.25,na.rm = T)-H) |
TRUE %in% (data>quantile(data,0.75,na.rm = T)+H))
outlier(data)
else
return(data)
}
Datos2.2_con_media<-Datos2.2
#Llamamos a la función de outliers para cada variable en la que hemos identificado un outlier
#Sería suficiente fijarnos en las gráficas donde hay outliers y no llamar a la función para todas las variables
Datos2.2_con_media[[6]]<-outlier(Datos2.2_con_media[[6]])
Datos2.2_con_media[[7]]<-outlier(Datos2.2_con_media[[7]])
Datos2.2_con_media[[8]]<-outlier(Datos2.2_con_media[[8]])
Datos2.2_con_media[[9]]<-outlier(Datos2.2_con_media[[9]])
Datos2.2_con_media[[10]]<-outlier(Datos2.2_con_media[[10]])
Datos2.2_con_media[[11]]<-outlier(Datos2.2_con_media[[11]])
Datos2.2_con_media[[12]]<-outlier(Datos2.2_con_media[[12]])
Datos2.2_con_media[[13]]<-outlier(Datos2.2_con_media[[13]])
Datos2.2_con_media[[14]]<-outlier(Datos2.2_con_media[[14]])
Datos2.2_con_media[[15]]<-outlier(Datos2.2_con_media[[15]])Comprobemos gráficamente el cambio de los outliers por la media.
par(mfrow = c(1, 2))
boxplot(Datos2.2[,6:15],main="Outliers antes de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
boxplot(Datos2.2_con_media[,6:15],main="Outliers después de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Estandarizamos.
par(mfrow = c(1, 2))
sca<-scale(Datos2.2[,6:15])
boxplot(sca,main="Outliers antes de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
sca_sin_outliers2<-scale(Datos2.2_con_media[,6:15])
boxplot(sca_sin_outliers,main="Outliers después de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Estudio 3
Eliminar los outliers
# Función para eliminar outliers de columnas específicas
eliminar_outliers_columnas <- function(datos, columnas) {
# Inicializar el dataframe sin outliers
datos_sin_outliers <- datos
for (columna in columnas) {
# Calcular los cuartiles y el IQR de la columna
Q1 <- quantile(datos[[columna]], 0.25)
Q3 <- quantile(datos[[columna]], 0.75)
IQR <- Q3 - Q1
# Definir los límites inferior y superior
limite_inferior <- Q1 - 1.5 * IQR
limite_superior <- Q3 + 1.5 * IQR
# Filtrar las filas que están dentro de los límites
datos_sin_outliers <- datos_sin_outliers[datos_sin_outliers[[columna]] >= limite_inferior & datos_sin_outliers[[columna]] <= limite_superior, ]
}
return(datos_sin_outliers)
}
columnas_a_usar <- c(6:15)
Datos3_sin_outliers <- eliminar_outliers_columnas(Datos3, columnas_a_usar)
#print(Datos3_sin_outliers)Comprobemos gráficamente la eliminación de los outliers.
par(mfrow = c(1, 2))
boxplot(Datos3[,6:15],main="Outliers antes de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
boxplot(Datos3_sin_outliers[,6:15],main="Outliers después de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Estandarizamos.
par(mfrow = c(1, 2))
sca<-scale(Datos3[,6:15])
boxplot(sca,main="Outliers antes de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
sca_sin_outliers<-scale(Datos3_sin_outliers[,6:15])
boxplot(sca_sin_outliers,main="Outliers después de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Cambiar los outliers por la media
outlier<-function(data,na.rm=T){
H<-1.5*IQR(data)
data[data<quantile(data,0.25,na.rm = T)-H]<-NA
data[data>quantile(data,0.75, na.rm = T)+H]<-NA
data[is.na(data)]<-mean(data, na.rm = T)
H<-1.5*IQR(data)
if (TRUE %in% (data<quantile(data,0.25,na.rm = T)-H) |
TRUE %in% (data>quantile(data,0.75,na.rm = T)+H))
outlier(data)
else
return(data)
}
Datos3_con_media<-Datos3
#Llamamos a la función de outliers para cada variable en la que hemos identificado un outlier
#Sería suficiente fijarnos en las gráficas donde hay outliers y no llamar a la función para todas las variables
Datos3_con_media[[6]]<-outlier(Datos3_con_media[[6]])
Datos3_con_media[[7]]<-outlier(Datos3_con_media[[7]])
Datos3_con_media[[8]]<-outlier(Datos3_con_media[[8]])
Datos3_con_media[[9]]<-outlier(Datos3_con_media[[9]])
Datos3_con_media[[10]]<-outlier(Datos3_con_media[[10]])
Datos3_con_media[[11]]<-outlier(Datos3_con_media[[11]])
Datos3_con_media[[12]]<-outlier(Datos3_con_media[[12]])
Datos3_con_media[[13]]<-outlier(Datos3_con_media[[13]])
Datos3_con_media[[14]]<-outlier(Datos3_con_media[[14]])
Datos3_con_media[[15]]<-outlier(Datos3_con_media[[15]])Comprobemos gráficamente el cambio de los outliers por la media.
par(mfrow = c(1, 2))
boxplot(Datos3[,6:15],main="Outliers antes de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
boxplot(Datos3_con_media[,6:15],main="Outliers después de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Estandarizamos.
par(mfrow = c(1, 2))
sca<-scale(Datos3[,6:15])
boxplot(sca,main="Outliers antes de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
sca_sin_outliers3<-scale(Datos3_con_media[,6:15])
boxplot(sca_sin_outliers,main="Outliers después de quitarlos",
xlab="All explanatory variables",
ylab="Values",
col=rainbow(10))
Análisis de Componentes Principales
Ahora, ya sin valores perdidos (los hemos eliminado), ni outliers (los hemos sustituido), hacemos ACP (Análisis de Componentes Principales)
Estudio de la correlación de las variables entre sí
Estudio 1
matriz_correlacion1<-cor(Datos1[6:15])
matriz_correlacion1## PesoFinal Tb.Sp Tb.N BMC BMD
## PesoFinal 1.0000000000 0.04899019 -0.076027758 0.88857725 0.67063992
## Tb.Sp 0.0489901901 1.00000000 -0.966650769 0.03772223 -0.07477822
## Tb.N -0.0760277581 -0.96665077 1.000000000 -0.07076835 0.06502617
## BMC 0.8885772453 0.03772223 -0.070768347 1.00000000 0.79200092
## BMD 0.6706399171 -0.07477822 0.065026167 0.79200092 1.00000000
## Longitud 0.7615998736 0.22750028 -0.274548124 0.82815136 0.58419435
## GPSurface 0.7479469259 0.16287707 -0.190656487 0.73317089 0.57053427
## GPVolume 0.3992803056 0.03884584 -0.008551150 0.20278668 0.16125363
## PesoInicial -0.0002011703 0.00913795 0.001712557 0.06261252 0.01841571
## PesoIntermedio 0.8614173499 0.11924157 -0.148865043 0.87152256 0.69088510
## Longitud GPSurface GPVolume PesoInicial
## PesoFinal 0.76159987 0.747946926 0.39928031 -0.0002011703
## Tb.Sp 0.22750028 0.162877069 0.03884584 0.0091379499
## Tb.N -0.27454812 -0.190656487 -0.00855115 0.0017125570
## BMC 0.82815136 0.733170887 0.20278668 0.0626125219
## BMD 0.58419435 0.570534274 0.16125363 0.0184157082
## Longitud 1.00000000 0.660639694 0.09404165 0.2078779916
## GPSurface 0.66063969 1.000000000 0.60361811 0.0034486662
## GPVolume 0.09404165 0.603618114 1.00000000 -0.1495347667
## PesoInicial 0.20787799 0.003448666 -0.14953477 1.0000000000
## PesoIntermedio 0.84301621 0.589832006 0.01870350 0.0512494835
## PesoIntermedio
## PesoFinal 0.86141735
## Tb.Sp 0.11924157
## Tb.N -0.14886504
## BMC 0.87152256
## BMD 0.69088510
## Longitud 0.84301621
## GPSurface 0.58983201
## GPVolume 0.01870350
## PesoInicial 0.05124948
## PesoIntermedio 1.00000000det(matriz_correlacion1)## [1] 1.311426e-05Estudio 2
matriz_correlacion2<-cor(Datos2.2[6:15])
matriz_correlacion2## PesoFinal Tb.Sp Tb.N BMC BMD
## PesoFinal 1.0000000 0.185923507 -0.21026340 0.928349174 0.81306413
## Tb.Sp 0.1859235 1.000000000 -0.96382803 -0.006849368 -0.07397760
## Tb.N -0.2102634 -0.963828027 1.00000000 -0.010983558 0.07778725
## BMC 0.9283492 -0.006849368 -0.01098356 1.000000000 0.89959506
## BMD 0.8130641 -0.073977602 0.07778725 0.899595064 1.00000000
## Longitud 0.8933615 0.081962798 -0.09378014 0.901471339 0.73605621
## GPSurface 0.7018824 0.058565382 -0.04145290 0.731245013 0.54036037
## GPVolume 0.1106489 -0.039730766 0.07318216 0.215404295 0.05951269
## PesoInicial 0.4002418 -0.028100641 0.05167504 0.442310869 0.37574638
## PesoIntermedio 0.9502286 0.248434899 -0.26003587 0.922008644 0.84160299
## Longitud GPSurface GPVolume PesoInicial PesoIntermedio
## PesoFinal 0.89336148 0.70188245 0.11064891 0.40024185 0.9502286
## Tb.Sp 0.08196280 0.05856538 -0.03973077 -0.02810064 0.2484349
## Tb.N -0.09378014 -0.04145290 0.07318216 0.05167504 -0.2600359
## BMC 0.90147134 0.73124501 0.21540430 0.44231087 0.9220086
## BMD 0.73605621 0.54036037 0.05951269 0.37574638 0.8416030
## Longitud 1.00000000 0.79971969 0.38438797 0.41533525 0.8214429
## GPSurface 0.79971969 1.00000000 0.64932125 0.10495674 0.6367105
## GPVolume 0.38438797 0.64932125 1.00000000 -0.09734002 0.0551226
## PesoInicial 0.41533525 0.10495674 -0.09734002 1.00000000 0.3858588
## PesoIntermedio 0.82144287 0.63671048 0.05512260 0.38585884 1.0000000det(matriz_correlacion2)## [1] 4.090502e-07Estudio 3
matriz_correlacion3<-cor(Datos3[6:15])
matriz_correlacion3## PesoFinal Tb.Sp Tb.N BMC BMD
## PesoFinal 1.00000000 -0.38603705 0.4068459 0.95790423 0.8590239
## Tb.Sp -0.38603705 1.00000000 -0.9723104 -0.50357610 -0.6405643
## Tb.N 0.40684592 -0.97231036 1.0000000 0.52230337 0.6476864
## BMC 0.95790423 -0.50357610 0.5223034 1.00000000 0.9322018
## BMD 0.85902394 -0.64056432 0.6476864 0.93220177 1.0000000
## Longitud 0.93550991 -0.35940363 0.3696785 0.94323641 0.8033721
## GPSurface 0.93119435 -0.32527309 0.3357800 0.90237031 0.8331614
## GPVolume 0.07754126 0.18163892 -0.1831888 0.05429013 0.0442623
## PesoInicial 0.13557640 -0.09253724 0.1228525 0.16003595 0.2106779
## PesoIntermedio 0.96651518 -0.35501921 0.3724031 0.94203518 0.8475858
## Longitud GPSurface GPVolume PesoInicial PesoIntermedio
## PesoFinal 0.93550991 0.9311944 0.07754126 0.13557640 0.96651518
## Tb.Sp -0.35940363 -0.3252731 0.18163892 -0.09253724 -0.35501921
## Tb.N 0.36967854 0.3357800 -0.18318880 0.12285251 0.37240306
## BMC 0.94323641 0.9023703 0.05429013 0.16003595 0.94203518
## BMD 0.80337207 0.8331614 0.04426230 0.21067786 0.84758578
## Longitud 1.00000000 0.8771301 0.02949129 0.14211526 0.92169706
## GPSurface 0.87713006 1.0000000 0.35807064 0.19751925 0.88068747
## GPVolume 0.02949129 0.3580706 1.00000000 0.25670125 -0.01902111
## PesoInicial 0.14211526 0.1975193 0.25670125 1.00000000 0.05779478
## PesoIntermedio 0.92169706 0.8806875 -0.01902111 0.05779478 1.00000000det(matriz_correlacion3)## [1] 1.816053e-08Test de Bartlett
Estudio 1
#dim(sca_sin_outliers1)
cortest.bartlett(cor(sca_sin_outliers1),n=56)## $chisq
## [1] 248.3985
##
## $p.value
## [1] 6.153321e-30
##
## $df
## [1] 45Estudio 2
#dim(sca_sin_outliers2)
cortest.bartlett(cor(sca_sin_outliers2),n=38)## $chisq
## [1] 84.19442
##
## $p.value
## [1] 0.0003579567
##
## $df
## [1] 45Estudio 3
#dim(sca_sin_outliers3)
cortest.bartlett(cor(sca_sin_outliers3),n=36)## $chisq
## [1] 442.8086
##
## $p.value
## [1] 8.5876e-67
##
## $df
## [1] 45Análisis de Componentes Principales
Estudio 1
PCA_1<-prcomp(Datos1_con_media[6:15], scale=T, center = T)
PCA_1$rotation## PC1 PC2 PC3 PC4 PC5
## PesoFinal -0.42237782 -0.01406992 0.18411363 0.306817447 -0.233732962
## Tb.Sp -0.17948406 -0.56845073 -0.35663205 -0.002435807 -0.006438170
## Tb.N 0.11605409 0.62190190 0.22929194 -0.020449368 -0.169758706
## BMC -0.43226109 0.19875483 0.01692992 0.111281614 0.009513443
## BMD -0.33313043 0.14521134 0.06551204 0.168566068 0.875819166
## Longitud -0.41128716 -0.02452259 -0.04627114 0.297477703 -0.273843241
## GPSurface -0.36371262 -0.21247245 0.34235604 -0.256417443 -0.159664452
## GPVolume -0.04739347 -0.25004961 0.65239142 -0.465139815 0.106184507
## PesoInicial -0.20799929 0.15643373 -0.46055366 -0.632033255 0.072553849
## PesoIntermedio -0.36212499 0.31076004 -0.15712924 -0.307533504 -0.179826720
## PC6 PC7 PC8 PC9 PC10
## PesoFinal 0.19230758 -0.070989722 -0.69511730 -0.21604990 0.25854548
## Tb.Sp -0.24792066 -0.243179428 0.12935726 0.21981729 0.57598751
## Tb.N -0.46797906 0.001188083 0.05570919 0.12103088 0.53028359
## BMC -0.04776174 -0.583724134 0.44868347 -0.44850427 -0.12408929
## BMD -0.07092137 0.138155692 -0.07161163 0.18693243 0.06311439
## Longitud 0.01301219 0.695216366 0.42168371 -0.04309907 0.03804733
## GPSurface -0.57986155 -0.014458407 -0.14486596 0.25849730 -0.43573252
## GPVolume 0.30788229 0.082702525 0.21554090 -0.16834823 0.32489087
## PesoInicial -0.13438586 0.248168413 -0.20638466 -0.43514547 0.06190064
## PesoIntermedio 0.47537135 -0.155047874 0.05171300 0.60653394 -0.01068467Vamos a calcular la desviación estándar de cada componente principal, la proporción de varianza explicada y la proporción de varianza explicada acumulada.
summary(PCA_1)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.9148 1.3710 1.1440 1.0408 0.80626 0.65430 0.61660
## Proportion of Variance 0.3666 0.1880 0.1309 0.1083 0.06501 0.04281 0.03802
## Cumulative Proportion 0.3666 0.5546 0.6855 0.7938 0.85881 0.90162 0.93964
## PC8 PC9 PC10
## Standard deviation 0.51334 0.45242 0.36792
## Proportion of Variance 0.02635 0.02047 0.01354
## Cumulative Proportion 0.96600 0.98646 1.00000La siguiente gráfica muestra la proporción de varianza explicada.
Explained_variance <- PCA_1$sdev^2 / sum(PCA_1$sdev^2)
ggplot(data = data.frame(Explained_variance, pc = 1:10),
aes(x = pc, y = Explained_variance, fill=Explained_variance )) +
geom_col(width = 0.3) + scale_y_continuous(limits = c(0,0.6)) + theme_bw() +
labs(x = "Principal component", y= "Proportion of variance")
La siguiente gráfica muestra la proporción de varianza explicada.
Cummulative_variance<-cumsum(Explained_variance)
ggplot( data = data.frame(Cummulative_variance, pc = 1:10),
aes(x = pc, y = Cummulative_variance ,fill=Cummulative_variance )) +
geom_col(width = 0.5) + scale_y_continuous(limits = c(0,1)) + theme_bw() +
labs(x = "Principal component",
y = "Proportion of cumulative variance")
Vemos cuántas componentes principales son adecuadas, usando la regla de Abdi: Las varianzas explicadas por los componentes principales se promedian, y se seleccionan aquellas cuya proporción de varianza explicada supera la media.
PCA_1$sdev^2## [1] 3.6664117 1.8796730 1.3087902 1.0832021 0.6500572 0.4281034 0.3801927
## [8] 0.2635212 0.2046870 0.1353615mean(PCA_1$sdev^2)## [1] 1Observamos que hay 4 por encima de la media, así que estudiamos el caso de 4 componentes principales. Esta gráfica muestra la proyección de las variables en dos dimensiones. Muestra el peso de la variable en la dirección de la componente principal.
# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca_var(PCA_1, axes = c(1, 2), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC1 vs PC2") + theme_bw()
# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca_var(PCA_1, axes = c(1, 3), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC1 vs PC3") + theme_bw()
# Gráfico 3: PC1 vs PC4
p3 <- fviz_pca_var(PCA_1, axes = c(1, 4), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC1 vs PC4") + theme_bw()
# Gráfico 4: PC2 vs PC3
p4 <- fviz_pca_var(PCA_1, axes = c(2, 3), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC2 vs PC3") + theme_bw()
# Gráfico 5: PC2 vs PC4
p5 <- fviz_pca_var(PCA_1, axes = c(2, 4), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC2 vs PC4") + theme_bw()
# Gráfico 6: PC3 vs PC4
p6 <- fviz_pca_var(PCA_1, axes = c(3, 4), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC3 vs PC4") + theme_bw()
grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)
Estos resultados gráficos muestran las observaciones con su contribución a la varianza. Así como identificar con colores aquellas observaciones que explican la mayor varianza de los componentes principales.
# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca_ind(PCA_1, col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, legend.title = "Contrib.var", title = "PC1 vs PC2") + theme_bw()
# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca_ind(PCA_1, axes = c(1, 3), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, legend.title = "Contrib.var", title = "PC1 vs PC3") + theme_bw()
# Gráfico 3: PC1 vs PC4
p3 <- fviz_pca_ind(PCA_1, axes = c(1, 4), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, legend.title = "Contrib.var", title = "PC1 vs PC4") + theme_bw()
# Gráfico 4: PC2 vs PC3
p4 <- fviz_pca_ind(PCA_1, axes = c(2, 3), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, legend.title = "Contrib.var", title = "PC2 vs PC3") + theme_bw()
# Gráfico 5: PC2 vs PC4
p5 <- fviz_pca_ind(PCA_1, axes = c(2, 4), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, legend.title = "Contrib.var", title = "PC2 vs PC4") + theme_bw()
# Gráfico 6: PC3 vs PC4
p6 <- fviz_pca_ind(PCA_1, axes = c(3, 4), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, legend.title = "Contrib.var", title = "PC3 vs PC4") + theme_bw()
grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)
Estos resultados gráficos muestran las variables y las observaciones con su contribución a la varianza. Identifica las posibles relaciones entre las contribuciones de los registros a las varianzas de los componentes y el peso de las variables en cada componente principal.
# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca(PCA_1, alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca(PCA_1, axes = c(1, 3), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
# Gráfico 3: PC1 vs PC4
p3 <- fviz_pca(PCA_1, axes = c(1, 4), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
# Gráfico 4: PC2 vs PC3
p4 <- fviz_pca(PCA_1, axes = c(2, 3), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
# Gráfico 5: PC2 vs PC4
p5 <- fviz_pca(PCA_1, axes = c(2, 4), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
# Gráfico 6: PC3 vs PC4
p6 <- fviz_pca(PCA_1, axes = c(3, 4), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)
Mostramos las coordenadas en el nuevo sistema de referencia
head(PCA_1$x)## PC1 PC2 PC3 PC4 PC5 PC6
## [1,] 0.2675163 -0.85141102 -0.155623192 -0.2840225 -1.2024258 0.27122398
## [2,] 2.8905478 -0.69771253 -0.117428533 0.9660797 -0.4560795 -0.84059860
## [3,] -2.3224893 -0.01339986 0.789150281 -0.6229904 0.5789967 0.29995198
## [4,] -1.3750855 1.37159909 0.008612358 -1.4599798 -1.1117157 0.93864317
## [5,] -3.0686746 0.89428864 -0.612825908 -0.5264444 0.1946402 0.04308362
## [6,] 1.4514801 -0.15087641 0.986361367 -0.1413722 -0.6230207 -0.30652335
## PC7 PC8 PC9 PC10
## [1,] -0.1082483 -0.06499945 -0.80442331 0.48813868
## [2,] -0.5972964 0.20512951 -0.35830257 -0.26738714
## [3,] -1.3578986 0.01939389 0.02074547 0.27015198
## [4,] 0.1181536 -0.47237701 0.55526997 0.31625801
## [5,] -0.6622062 0.13964267 -0.59058134 -0.07965825
## [6,] -0.2994882 0.39283662 0.08906888 -0.10630711Estudio 2
PCA_2<-prcomp(Datos2.2_con_media[6:15], scale=T, center = T)
PCA_2$rotation## PC1 PC2 PC3 PC4 PC5
## PesoFinal -0.25992492 0.2980351 -0.34757661 -0.47859323 0.17107702
## Tb.Sp -0.23710824 -0.3098604 0.28815744 -0.40296686 0.43989353
## Tb.N 0.06134366 0.4273067 -0.42280589 0.35261251 -0.01130212
## BMC 0.33185093 0.1541500 -0.08552962 0.11718384 0.78389868
## BMD -0.50447365 -0.1280594 -0.19453870 0.33384119 -0.01191529
## Longitud 0.33859316 0.2399386 -0.24789866 -0.53977919 -0.24429972
## GPSurface -0.14483379 -0.4270193 -0.47112369 -0.22240881 -0.12003660
## GPVolume 0.12791088 -0.3892311 -0.52250284 0.11078784 0.24968223
## PesoInicial -0.43493454 0.1285923 -0.10430309 0.01505287 0.04648235
## PesoIntermedio -0.40827143 0.4294838 0.07021918 -0.07563641 0.15448793
## PC6 PC7 PC8 PC9 PC10
## PesoFinal 0.22304724 0.02845598 -0.395331819 0.41670553 0.289819325
## Tb.Sp 0.08193385 0.44118628 0.451740845 0.04821205 -0.033018387
## Tb.N 0.14110777 0.44296051 0.483238167 -0.02806575 0.244919554
## BMC -0.13025436 -0.44156481 0.124754826 0.05528321 0.001598339
## BMD 0.17456712 -0.25400188 0.213287489 0.46060872 -0.476473841
## Longitud -0.17866716 -0.04959764 0.291311656 0.06432760 -0.540281481
## GPSurface 0.06163855 -0.41281942 0.304729359 -0.34523956 0.352953462
## GPVolume -0.02092199 0.39980001 -0.400720968 -0.20776327 -0.347120891
## PesoInicial -0.87666584 0.07513788 -0.009934938 -0.04509557 0.068619429
## PesoIntermedio 0.26736399 -0.10231783 -0.081670411 -0.66294812 -0.293627855Vamos a calcular la desviación estándar de cada componente principal, la proporción de varianza explicada y la proporción de varianza explicada acumulada.
summary(PCA_2)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.5052 1.3747 1.2400 1.1146 0.96014 0.81996 0.68561
## Proportion of Variance 0.2266 0.1890 0.1538 0.1242 0.09219 0.06723 0.04701
## Cumulative Proportion 0.2266 0.4156 0.5693 0.6935 0.78574 0.85297 0.89998
## PC8 PC9 PC10
## Standard deviation 0.65022 0.61002 0.45311
## Proportion of Variance 0.04228 0.03721 0.02053
## Cumulative Proportion 0.94226 0.97947 1.00000La siguiente gráfica muestra la proporción de varianza explicada.
Explained_variance <- PCA_2$sdev^2 / sum(PCA_2$sdev^2)
ggplot(data = data.frame(Explained_variance, pc = 1:10),
aes(x = pc, y = Explained_variance, fill=Explained_variance )) +
geom_col(width = 0.3) + scale_y_continuous(limits = c(0,0.6)) + theme_bw() +
labs(x = "Principal component", y= "Proportion of variance")
La siguiente gráfica muestra la proporción de varianza explicada.
Cummulative_variance<-cumsum(Explained_variance)
ggplot( data = data.frame(Cummulative_variance, pc = 1:10),
aes(x = pc, y = Cummulative_variance ,fill=Cummulative_variance )) +
geom_col(width = 0.5) + scale_y_continuous(limits = c(0,1)) + theme_bw() +
labs(x = "Principal component",
y = "Proportion of cumulative variance")
Vemos cuántas componentes principales son adecuadas, usando la regla de Abdi: Las varianzas explicadas por los componentes principales se promedian, y se seleccionan aquellas cuya proporción de varianza explicada supera la media.
PCA_2$sdev^2## [1] 2.2657583 1.8898238 1.5376064 1.2423200 0.9218780 0.6723380 0.4700644
## [8] 0.4227817 0.3721201 0.2053092mean(PCA_2$sdev^2)## [1] 1Observamos que hay 4 por encima de la media, así que estudiamos el caso de 4 componentes principales. Esta gráfica muestra la proyección de las variables en dos dimensiones. Muestra el peso de la variable en la dirección de la componente principal.
# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca_var(PCA_2, axes = c(1, 2), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC1 vs PC2") + theme_bw()
# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca_var(PCA_2, axes = c(1, 3), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC1 vs PC3") + theme_bw()
# Gráfico 3: PC1 vs PC4
p3 <- fviz_pca_var(PCA_2, axes = c(1, 4), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC1 vs PC4") + theme_bw()
# Gráfico 4: PC2 vs PC3
p4 <- fviz_pca_var(PCA_2, axes = c(2, 3), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC2 vs PC3") + theme_bw()
# Gráfico 5: PC2 vs PC4
p5 <- fviz_pca_var(PCA_2, axes = c(2, 4), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC2 vs PC4") + theme_bw()
# Gráfico 6: PC3 vs PC4
p6 <- fviz_pca_var(PCA_2, axes = c(3, 4), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC3 vs PC4") + theme_bw()
grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)
Estos resultados gráficos muestran las observaciones con su contribución a la varianza. Así como identificar con colores aquellas observaciones que explican la mayor varianza de los componentes principales.
# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca_ind(PCA_2, col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, legend.title = "Contrib.var", title = "PC1 vs PC2") + theme_bw()
# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca_ind(PCA_2, axes = c(1, 3), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, legend.title = "Contrib.var", title = "PC1 vs PC3") + theme_bw()
# Gráfico 3: PC1 vs PC4
p3 <- fviz_pca_ind(PCA_2, axes = c(1, 4), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, legend.title = "Contrib.var", title = "PC1 vs PC4") + theme_bw()
# Gráfico 4: PC2 vs PC3
p4 <- fviz_pca_ind(PCA_2, axes = c(2, 3), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, legend.title = "Contrib.var", title = "PC2 vs PC3") + theme_bw()
# Gráfico 5: PC2 vs PC4
p5 <- fviz_pca_ind(PCA_2, axes = c(2, 4), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, legend.title = "Contrib.var", title = "PC2 vs PC4") + theme_bw()
# Gráfico 6: PC3 vs PC4
p6 <- fviz_pca_ind(PCA_2, axes = c(3, 4), col.ind = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, legend.title = "Contrib.var", title = "PC3 vs PC4") + theme_bw()
grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)
Estos resultados gráficos muestran las variables y las observaciones con su contribución a la varianza. Identifica las posibles relaciones entre las contribuciones de los registros a las varianzas de los componentes y el peso de las variables en cada componente principal.
# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca(PCA_2, alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca(PCA_2, axes = c(1, 3), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
# Gráfico 3: PC1 vs PC4
p3 <- fviz_pca(PCA_2, axes = c(1, 4), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
# Gráfico 4: PC2 vs PC3
p4 <- fviz_pca(PCA_2, axes = c(2, 3), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
# Gráfico 5: PC2 vs PC4
p5 <- fviz_pca(PCA_2, axes = c(2, 4), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
# Gráfico 6: PC3 vs PC4
p6 <- fviz_pca(PCA_2, axes = c(3, 4), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 3)
Mostramos las coordenadas en el nuevo sistema de referencia
head(PCA_2$x)## PC1 PC2 PC3 PC4 PC5 PC6
## [1,] -0.5565609 1.0555072 -0.5221437 -0.14240038 -0.03484642 0.9084108
## [2,] 0.9333002 -0.9384008 2.7251929 -0.07547032 -0.89002717 0.4401783
## [3,] 2.2820594 0.6146465 2.1675107 -0.56366156 -0.12542127 -0.2913732
## [4,] -1.0898128 1.4132247 -1.9715946 0.53085606 0.96370299 0.6789528
## [5,] -0.5830088 1.2057056 1.0839525 -0.89863163 0.49375772 0.1674253
## [6,] -0.5097828 -0.9369064 0.5602292 -2.34375552 -1.54888116 0.4859421
## PC7 PC8 PC9 PC10
## [1,] -0.46745970 -0.31353621 -0.5008955 0.0313710
## [2,] -0.30552079 -0.06875839 -1.1897007 0.9964603
## [3,] -0.63626742 0.34002247 0.1058022 -0.7520804
## [4,] -0.04304059 1.12403097 -0.8339104 0.6856212
## [5,] 0.33700562 -0.01534292 -0.1676559 0.2327255
## [6,] 1.56762171 0.76554558 -0.3195342 -0.3117803Estudio 3
PCA_3<-prcomp(Datos3_con_media[6:15], scale=T, center = T)
PCA_3$rotation## PC1 PC2 PC3 PC4 PC5
## PesoFinal -0.39673358 -0.15725522 -0.05974109 -0.02209423 0.223708700
## Tb.Sp 0.25415886 -0.54303826 -0.22917115 -0.10898013 0.269494012
## Tb.N -0.25208263 0.54499545 0.19743935 0.24929268 -0.003595032
## BMC -0.40920080 -0.06996250 -0.02757952 0.02344307 0.195991225
## BMD -0.39659089 0.05588512 0.04864212 0.12142489 -0.039844026
## Longitud -0.38263075 -0.16813439 -0.12425103 -0.07929067 0.365364278
## GPSurface -0.37830852 -0.26824077 0.17112292 0.01737851 0.046277032
## GPVolume -0.02017475 -0.42080586 0.79776660 -0.02411437 -0.296473049
## PesoInicial 0.14501040 -0.21612759 -0.03968514 0.94702119 0.123818798
## PesoIntermedio -0.28066312 -0.22656417 -0.46805425 0.07858178 -0.773598333
## PC6 PC7 PC8 PC9 PC10
## PesoFinal 0.228149851 -0.235384793 -0.49970271 -0.522602781 -0.365370643
## Tb.Sp -0.439670214 0.437303981 -0.33921780 0.008999114 0.061605549
## Tb.N -0.067634099 0.624420509 -0.37609358 0.007990931 0.048027822
## BMC -0.147175330 0.001963112 0.28414183 -0.418533383 0.714399034
## BMD -0.776490179 -0.253618142 0.13139464 0.143887971 -0.340099181
## Longitud 0.284816556 0.417927267 0.49451461 0.221286421 -0.347580355
## GPSurface 0.175315256 -0.213181226 -0.36457661 0.664063133 0.321131282
## GPVolume 0.004249656 0.192082291 0.11074756 -0.196659068 -0.098755359
## PesoInicial 0.090358748 -0.059573312 0.07651613 -0.016323749 -0.025900613
## PesoIntermedio 0.073669048 0.196413387 -0.01415723 -0.044483427 0.002842317Vamos a calcular la desviación estándar de cada componente principal, la proporción de varianza explicada y la proporción de varianza explicada acumulada.
summary(PCA_3)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.3934 1.2892 1.0217 0.94026 0.61297 0.3377 0.31351
## Proportion of Variance 0.5728 0.1662 0.1044 0.08841 0.03757 0.0114 0.00983
## Cumulative Proportion 0.5728 0.7390 0.8434 0.93181 0.96938 0.9808 0.99061
## PC8 PC9 PC10
## Standard deviation 0.23525 0.16921 0.09942
## Proportion of Variance 0.00553 0.00286 0.00099
## Cumulative Proportion 0.99615 0.99901 1.00000La siguiente gráfica muestra la proporción de varianza explicada.
Explained_variance <- PCA_3$sdev^2 / sum(PCA_3$sdev^2)
ggplot(data = data.frame(Explained_variance, pc = 1:10),
aes(x = pc, y = Explained_variance, fill=Explained_variance )) +
geom_col(width = 0.3) + scale_y_continuous(limits = c(0,0.6)) + theme_bw() +
labs(x = "Principal component", y= "Proportion of variance")
La siguiente gráfica muestra la proporción de varianza explicada.
Cummulative_variance<-cumsum(Explained_variance)
ggplot( data = data.frame(Cummulative_variance, pc = 1:10),
aes(x = pc, y = Cummulative_variance ,fill=Cummulative_variance )) +
geom_col(width = 0.5) + scale_y_continuous(limits = c(0,1)) + theme_bw() +
labs(x = "Principal component",
y = "Proportion of cumulative variance")
Vemos cuántas componentes principales son adecuadas, usando la regla de Abdi: Las varianzas explicadas por los componentes principales se promedian, y se seleccionan aquellas cuya proporción de varianza explicada supera la media.
PCA_3$sdev^2## [1] 5.728162672 1.662065418 1.043793886 0.884082735 0.375734284 0.114015776
## [7] 0.098286304 0.055342812 0.028632712 0.009883402mean(PCA_3$sdev^2)## [1] 1Observamos que hay 3 por encima de la media, así que estudiamos el caso de 3 componentes principales. Esta gráfica muestra la proyección de las variables en dos dimensiones. Muestra el peso de la variable en la dirección de la componente principal.
# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca_var(PCA_3, axes = c(1, 2), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC1 vs PC2") + theme_bw()
# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca_var(PCA_3, axes = c(1, 3), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC1 vs PC3") + theme_bw()
# Gráfico 3: PC2 vs PC3
p3 <- fviz_pca_var(PCA_3, axes = c(2, 3), repel = TRUE, col.var = "cos2",
legend.title = "Distance", title = "PC2 vs PC3") + theme_bw()
grid.arrange(p1, p2, p3, nrow = 2)
Estos resultados gráficos muestran las observaciones con su contribución a la varianza. Así como identificar con colores aquellas observaciones que explican la mayor varianza de los componentes principales.
# Gráfico 1: PC1 vs PC2
p1<-fviz_pca_ind(PCA_3,col.ind = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel=TRUE,legend.title="Contrib.var", title="Records")+theme_bw()
# Gráfico 2: PC1 vs PC3
p2<-fviz_pca_ind(PCA_3,axes=c(1,3),col.ind = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel=TRUE,legend.title="Contrib.var", title="Records")+theme_bw()
# Gráfico 1: PC2 vs PC3
p3<-fviz_pca_ind(PCA_3,axes=c(2,3),col.ind = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel=TRUE,legend.title="Contrib.var", title="Records")+theme_bw()
grid.arrange(p1, p2, p3, nrow = 2)
Estos resultados gráficos muestran las variables y las observaciones con su contribución a la varianza. Identifica las posibles relaciones entre las contribuciones de los registros a las varianzas de los componentes y el peso de las variables en cada componente principal.
# Gráfico 1: PC1 vs PC2
p1 <- fviz_pca(PCA_3, alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
# Gráfico 2: PC1 vs PC3
p2 <- fviz_pca(PCA_3, axes = c(1, 3), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
# Gráfico 3: PC2 vs PC3
p3 <- fviz_pca(PCA_3, axes = c(2, 3), alpha.ind = "contrib", col.var = "cos2", col.ind = "seagreen",
gradient.cols = c("#FDF50E", "#FD960E", "#FD1E0E"), repel = TRUE, legend.title = "Distancia") + theme_bw()
grid.arrange(p1, p2, p3, nrow = 2)
Mostramos las coordenadas en el nuevo sistema de referencia
head(PCA_3$x)## PC1 PC2 PC3 PC4 PC5 PC6
## [1,] -1.5233677 -0.3581962 0.9210322 0.41030042 0.12089330 0.02626068
## [2,] 3.1185773 0.6346817 -0.5727801 0.03208267 -0.96206638 -0.33121675
## [3,] 1.0295563 0.4364019 1.4661299 0.54437498 0.01583101 0.29865500
## [4,] 0.5514543 0.9841047 0.3842833 -0.81213671 -0.66415803 -0.25097263
## [5,] 1.2337493 0.9170182 0.3097037 0.75243837 0.16870026 0.39776375
## [6,] 0.4911084 1.0929372 -0.4006754 0.16191980 -0.08825771 0.25773892
## PC7 PC8 PC9 PC10
## [1,] 0.06277742 0.09463181 0.18848194 -0.15272820
## [2,] 0.20938983 0.22941818 -0.25554403 -0.07244920
## [3,] -0.22116901 0.02723566 0.13598678 0.05547987
## [4,] 0.16641356 0.12817490 0.02373183 0.02365107
## [5,] -0.05732605 0.02661690 0.26384796 -0.09626547
## [6,] -0.01416743 -0.21702662 -0.02839555 0.05045322#print(PCA_3$x)
#print(PCA_3)Estudio de la normalidad
Previo a la construcción de métodos de clasificación se debe analizar la normalidad multivariante de los datos con el test propuesto en clase de prácticas de análisis discriminante.
Estudio 1
PCA_1scores<-as.data.frame(PCA_1$x[,1:4])
outliers1 <- mvn(data = PCA_1scores, mvnTest = "hz", multivariateOutlierMethod = "quan")
Estudio 2
PCA_2scores<-as.data.frame(PCA_2$x[,1:4])
outliers2 <- mvn(data = PCA_2scores, mvnTest = "hz", multivariateOutlierMethod = "quan")
Estudio 3
PCA_3scores<-as.data.frame(PCA_3$x[,1:3])
outliers3 <- mvn(data = PCA_3scores, mvnTest = "hz", multivariateOutlierMethod = "quan")
Royston multivariate normality test
Estudio 1
royston_test1 <- mvn(data = PCA_1scores, mvnTest = "royston", multivariatePlot = "qq")
royston_test1$multivariateNormality## Test H p value MVN
## 1 Royston 17.61511 0.001467199 NOEstudio 2
royston_test2 <- mvn(data = PCA_2scores, mvnTest = "royston", multivariatePlot = "qq")
royston_test2$multivariateNormality## Test H p value MVN
## 1 Royston 2.532671 0.6387953 YESEstudio 3
royston_test3 <- mvn(data = PCA_3scores, mvnTest = "royston", multivariatePlot = "qq")
royston_test3$multivariateNormality## Test H p value MVN
## 1 Royston 10.02739 0.01833478 NOHenze-Zirkler multivariate normality test
Estudio 1
hz_test1 <- mvn(data = PCA_1scores, mvnTest = "hz")
hz_test1$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 0.929499 0.07517825 YESEstudio 2
hz_test2 <- mvn(data = PCA_2scores, mvnTest = "hz")
hz_test2$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 0.811583 0.2443798 YESEstudio 3
hz_test3 <- mvn(data = PCA_3scores, mvnTest = "hz")
hz_test3$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 1.246543 0.000942336 NOLas variables respuesta (dependientes) serían: Tb.Sp, Tb.N, BMC, BMD, Longitud, GPSurface, GPVolume.
Las variables independientes (predictoras) serían: Grupo, Dosis, DiasDosis,Peso, PesoInicial, Peso Intermedio, DiasEspera.
Análisis clúster
Reescalamos las componentes principales.
PCA_1scores_scale<-PCA_1scores
PCA_1scores_scale<-as.data.frame(scale(PCA_1scores_scale))
PCA_2scores_scale<-PCA_2scores
PCA_2scores_scale<-as.data.frame(scale(PCA_2scores_scale))
PCA_3scores_scale<-PCA_3scores
PCA_3scores_scale<-as.data.frame(scale(PCA_3scores_scale))Hacemos la matriz de distancia para cada estudio.
Estudio 1
distance<- get_dist(PCA_1scores_scale)
fviz_dist(distance, gradient = list(low ="yellow", mid = "white", high = "brown"))
Estudio 2
distance<- get_dist(PCA_2scores_scale)
fviz_dist(distance, gradient = list(low ="yellow", mid = "white", high = "brown"))
Estudio 3
distance<- get_dist(PCA_3scores_scale)
fviz_dist(distance, gradient = list(low ="yellow", mid = "white", high = "brown"))
Clustering jerárquico - método de Ward
Hacemos los dendogramas.
Estudio 1
dendrogram <- hclust(dist(PCA_1scores_scale, method = 'euclidean'), method = 'ward.D')
ggdendrogram(dendrogram, rotate = FALSE, labels = FALSE, theme_dendro = TRUE) +
labs(title = "Dendrograma")
Estudio 2
dendrogram <- hclust(dist(PCA_2scores_scale, method = 'euclidean'), method = 'ward.D')
ggdendrogram(dendrogram, rotate = FALSE, labels = FALSE, theme_dendro = TRUE) +
labs(title = "Dendrograma")
Estudio 3
dendrogram <- hclust(dist(PCA_3scores_scale, method = 'euclidean'), method = 'ward.D')
ggdendrogram(dendrogram, rotate = FALSE, labels = FALSE, theme_dendro = TRUE) +
labs(title = "Dendrograma")
Clustering no jerárquico - algoritmo K-means
Estudio 1
set.seed(123)
k2_1 <- kmeans(PCA_1scores_scale, centers = 2, nstart = 25)
k3_1 <- kmeans(PCA_1scores_scale, centers = 3, nstart = 25)
k4_1 <- kmeans(PCA_1scores_scale, centers = 4, nstart = 25)
k5_1 <- kmeans(PCA_1scores_scale, centers = 5, nstart = 25)
k6_1 <- kmeans(PCA_1scores_scale, centers = 6, nstart = 25)
k7_1 <- kmeans(PCA_1scores_scale, centers = 7, nstart = 25)
k9_1 <- kmeans(PCA_1scores_scale, centers = 9, nstart = 25)
# Plots to compare
fviz_cluster(k2_1, data = PCA_1scores_scale) + ggtitle("k = 2") 
fviz_cluster(k3_1, data = PCA_1scores_scale) + ggtitle("k = 3") 
fviz_cluster(k4_1, data = PCA_1scores_scale) + ggtitle("k = 4") 
fviz_cluster(k5_1, data = PCA_1scores_scale) + ggtitle("k = 5")
fviz_cluster(k6_1, data = PCA_1scores_scale) + ggtitle("k = 6") 
fviz_cluster(k7_1, data = PCA_1scores_scale) + ggtitle("k = 7")
fviz_cluster(k9_1, data = PCA_1scores_scale) + ggtitle("k = 9") 
Apliquemos el método de Elbow, de Silhouette y el del estadístico de brecha para ver el número óptimo de clústers.
Método de Elbow
set.seed(123)
fviz_nbclust(PCA_1scores_scale, kmeans, method = "wss")
Método de Silhouette
set.seed(123)
fviz_nbclust(PCA_1scores_scale, kmeans, method = "silhouette")
Método estadístico de brecha (GAP)
set.seed(123)
gap_stat <- clusGap(PCA_1scores_scale, FUN = kmeans, nstart = 25,K.max = 10, B = 50)
fviz_gap_stat(gap_stat)
Estudio 2
set.seed(123)
k2_2 <- kmeans(PCA_2scores_scale, centers = 2, nstart = 25)
k3_2 <- kmeans(PCA_2scores_scale, centers = 3, nstart = 25)
k4_2 <- kmeans(PCA_2scores_scale, centers = 4, nstart = 25)
k5_2 <- kmeans(PCA_2scores_scale, centers = 5, nstart = 25)
# Plots to compare
fviz_cluster(k2_2, data = PCA_2scores_scale) + ggtitle("k = 2") 
fviz_cluster(k3_2, data = PCA_2scores_scale) + ggtitle("k = 3") 
fviz_cluster(k4_2, data = PCA_2scores_scale) + ggtitle("k = 4") 
fviz_cluster(k5_2, data = PCA_2scores_scale) + ggtitle("k = 5")
Apliquemos el método de Elbow, de Silhouette y el del estadístico de brecha para ver el número óptimo de clústers.
Método de Elbow
set.seed(123)
fviz_nbclust(PCA_2scores_scale, kmeans, method = "wss")
Método de Silhouette
set.seed(123)
fviz_nbclust(PCA_2scores_scale, kmeans, method = "silhouette")
Método estadístico de brecha (GAP)
set.seed(123)
gap_stat <- clusGap(PCA_2scores_scale, FUN = kmeans, nstart = 25,K.max = 10, B = 50)
fviz_gap_stat(gap_stat)
Estudio 3 ( sin escalar )
set.seed(123)
k2_3 <- kmeans(PCA_3scores, centers = 2, nstart = 25)
k3_3 <- kmeans(PCA_3scores, centers = 3, nstart = 25)
k4_3 <- kmeans(PCA_3scores, centers = 4, nstart = 25)
k5_3 <- kmeans(PCA_3scores, centers = 5, nstart = 25)
# Plots to compare
fviz_cluster(k2_3, data = PCA_3scores) + ggtitle("k = 2")
fviz_cluster(k3_3, data = PCA_3scores) + ggtitle("k = 3")
fviz_cluster(k4_3, data = PCA_3scores) + ggtitle("k = 4")
fviz_cluster(k5_3, data = PCA_3scores) + ggtitle("k = 5")
Apliquemos el método de Elbow, de Silhouette y el del estadístico de brecha para ver el número óptimo de clústers.
Método de Elbow
set.seed(123)
fviz_nbclust(PCA_3scores_scale, kmeans, method = "wss")
Método de Silhouette
set.seed(123)
fviz_nbclust(PCA_3scores_scale, kmeans, method = "silhouette")
Método estadístico de brecha (GAP)
set.seed(123)
gap_stat <- clusGap(PCA_3scores_scale, FUN = kmeans, nstart = 25,K.max = 10, B = 50)
fviz_gap_stat(gap_stat)
Cálculo de vectores de media
Veamos el número óptimo de clusters. Calculamos vectores de medias de 2, 3 y 4 clusters. #### Estudio 1 ##### Vectores de medias con 2 cluster.
Datos1_con_media$Cluster2 <- k2_1$cluster
Datos1_con_media$Cluster3 <- k3_1$cluster
Datos1_con_media$Cluster4 <- k4_1$cluster
colMeans(subset(Datos1_con_media[3:17], Cluster2 == 1))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.0000000 21.1000000 14.0000000 288.5552830 0.3245529
## Tb.N BMC BMD Longitud GPSurface
## 2.9881437 0.2632654 142.2262667 38.7577212 39.1877675
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 4.7023645 123.1900000 164.6585354 14.0000000 1.0000000colMeans(subset(Datos1_con_media[3:17], Cluster2 == 2))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.0000000 36.8846154 14.0000000 304.2852685 0.2933026
## Tb.N BMC BMD Longitud GPSurface
## 3.2657038 0.2682001 143.4288177 38.9815385 41.2414331
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.9588533 116.5961538 165.1677350 14.0000000 2.0000000Vectores de medias con 3 cluster.
colMeans(subset(Datos1_con_media[3:18], Cluster3 == 1))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.0000000 25.7600000 14.0000000 277.7183396 0.3170226
## Tb.N BMC BMD Longitud GPSurface
## 2.9973644 0.2509190 138.9594800 38.2376655 38.5942069
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.0396320 118.5800000 158.3110480 14.0000000 1.2000000
## Cluster3
## 1.0000000colMeans(subset(Datos1_con_media[3:18], Cluster3 == 2))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.0000000 36.4615385 14.0000000 303.2198839 0.2967179
## Tb.N BMC BMD Longitud GPSurface
## 3.3144923 0.2766495 145.1298718 39.0961538 41.5172458
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.7030430 122.4230769 171.1128739 14.0000000 1.7307692
## Cluster3
## 2.0000000colMeans(subset(Datos1_con_media[3:18], Cluster3 == 3))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.0000000 0.0000000 14.0000000 348.2800000 0.3444445
## Tb.N BMC BMD Longitud GPSurface
## 2.6883400 0.2810602 149.7147188 40.7620000 40.7213439
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 4.3462400 115.9400000 165.4812500 14.0000000 1.4000000
## Cluster3
## 3.0000000colMeans(subset(Datos1_con_media[3:18], Dosis == 80))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.0000000 80.0000000 14.0000000 285.6750000 0.2910465
## Tb.N BMC BMD Longitud GPSurface
## 3.3513250 0.2647625 142.1288750 38.7700000 41.3159225
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.6107819 119.7375000 162.1545139 14.0000000 1.6250000
## Cluster3
## 1.7500000Vectores de medias con 4 cluster.
colMeans(subset(Datos1_con_media[3:19], Cluster4 == 1))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.0000000 21.7000000 14.0000000 275.2250000 0.3052372
## Tb.N BMC BMD Longitud GPSurface
## 3.1449455 0.2545940 139.7934000 38.3180000 37.8528443
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 4.6114900 120.5550000 161.5147475 14.0000000 1.0000000
## Cluster3 Cluster4
## 1.1000000 1.0000000colMeans(subset(Datos1_con_media[3:19], Cluster4 == 2))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.0000000 39.4210526 14.0000000 295.2587885 0.2854012
## Tb.N BMC BMD Longitud GPSurface
## 3.4086947 0.2624125 141.9950877 38.6578947 40.8136171
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 6.0205033 115.6631579 161.9559211 14.0000000 2.0000000
## Cluster3 Cluster4
## 1.7368421 2.0000000colMeans(subset(Datos1_con_media[3:19], Cluster4 == 3))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.0000000 0.0000000 14.0000000 346.6250000 0.3304307
## Tb.N BMC BMD Longitud GPSurface
## 2.7367000 0.2849361 149.0038985 40.7675000 40.0302841
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 4.3432750 109.7250000 165.4812500 14.0000000 1.5000000
## Cluster3 Cluster4
## 3.0000000 3.0000000colMeans(subset(Datos1_con_media[3:19], Cluster4 == 4))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.000000 31.461538 14.000000 312.858345 0.347182
## Tb.N BMC BMD Longitud GPSurface
## 2.764746 0.281054 146.626692 39.409357 42.713349
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.539127 129.200000 174.210363 14.000000 1.384615
## Cluster3 Cluster4
## 1.923077 4.000000Estudio 2
Vectores de medias con 2 cluster.
Datos2.2_con_media$Cluster2 <- k2_2$cluster
Datos2.2_con_media$Cluster3 <- k3_2$cluster
Datos2.2_con_media$Cluster4 <- k4_2$cluster
colMeans(subset(Datos2.2_con_media[3:17], Cluster2 == 1))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 176.8421053 21.0000000 351.7894737 0.3148076
## Tb.N BMC BMD Longitud GPSurface
## 3.1184895 0.3245605 154.3637895 41.7489522 43.2144266
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.4325263 115.9684211 164.3157895 28.0000000 1.0000000colMeans(subset(Datos2.2_con_media[3:17], Cluster2 == 2))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 58.9473684 21.0000000 331.9638596 0.3190193
## Tb.N BMC BMD Longitud GPSurface
## 2.8985947 0.3242006 164.3222632 41.7358057 43.8724149
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.8739789 116.2157895 154.9175439 28.0000000 2.0000000Vectores de medias con 3 cluster.
colMeans(subset(Datos2.2_con_media[3:18], Cluster3 == 1))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 175.0000000 21.0000000 356.4312500 0.3167109
## Tb.N BMC BMD Longitud GPSurface
## 3.1983813 0.3244072 156.0120000 41.7480842 43.2287417
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.4088875 117.3312500 169.1937500 28.0000000 1.0000000
## Cluster3
## 1.0000000colMeans(subset(Datos2.2_con_media[3:18], Cluster3 == 2))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 8.8888889 21.0000000 341.0237037 0.3253313
## Tb.N BMC BMD Longitud GPSurface
## 2.7106444 0.3234977 176.6955556 41.7344344 45.6307988
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.9622333 120.6555556 161.0037037 28.0000000 2.0000000
## Cluster3
## 2.0000000colMeans(subset(Datos2.2_con_media[3:18], Cluster3 == 3))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.000000 123.076923 21.000000 324.553846 0.311335
## Tb.N BMC BMD Longitud GPSurface
## 2.981131 0.324959 151.429462 41.740857 42.485610
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.740100 111.407692 146.869231 28.000000 1.769231
## Cluster3
## 3.000000Vectores de medias con 4 cluster.
colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 1))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 8.8888889 21.0000000 341.0237037 0.3253313
## Tb.N BMC BMD Longitud GPSurface
## 2.7106444 0.3234977 176.6955556 41.7344344 45.6307988
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.9622333 120.6555556 161.0037037 28.0000000 2.0000000
## Cluster3 Cluster4
## 2.0000000 1.0000000colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 2))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 160.0000000 21.0000000 354.6500000 0.3189563
## Tb.N BMC BMD Longitud GPSurface
## 3.1786250 0.3242581 157.4050000 41.7466503 42.7150195
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.0384500 117.8583333 175.7666667 28.0000000 1.0000000
## Cluster3 Cluster4
## 1.0000000 2.0000000colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 3))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 70.0000000 21.0000000 321.5750000 0.3133542
## Tb.N BMC BMD Longitud GPSurface
## 2.8665125 0.3247514 151.8255000 41.7371055 41.3465789
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.1177875 111.9750000 151.6625000 28.0000000 1.8750000
## Cluster3 Cluster4
## 3.0000000 3.0000000colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 4))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 213.3333333 21.0000000 343.7444444 0.3089356
## Tb.N BMC BMD Longitud GPSurface
## 3.2059111 0.3250970 151.2567778 41.7493159 44.5133260
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 6.6399778 112.8333333 143.7666667 28.0000000 1.3333333
## Cluster3 Cluster4
## 2.1111111 4.0000000colMeans(subset(Datos2.2_con_media[3:19], Dosis == 80))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 80.0000000 21.0000000 339.7090909 0.3234777
## Tb.N BMC BMD Longitud GPSurface
## 3.1073455 0.3247939 153.7318182 41.7446471 44.1898289
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.8196273 113.6636364 158.6727273 28.0000000 1.3636364
## Cluster3 Cluster4
## 2.0000000 2.8181818Estudio 3
Vectores de medias con 2 cluster.
Datos3_con_media$Cluster2 <- k2_3$cluster
Datos3_con_media$Cluster3 <- k3_3$cluster
Datos3_con_media$Cluster4 <- k4_3$cluster
Datos3_con_media$Cluster5 <- k5_3$cluster
colMeans(subset(Datos3_con_media[3:17], Cluster2 == 1))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 3.0000000 519.2307692 14.0000000 211.9230769 0.3209199
## Tb.N BMC BMD Longitud GPSurface
## 3.1546231 0.1793077 118.0496923 33.9794231 33.9432231
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.6273692 71.7447653 99.1657484 14.0000000 1.0000000colMeans(subset(Datos3_con_media[3:17], Cluster2 == 2))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 3.00000 0.00000 14.00000 306.80000 0.27303
## Tb.N BMC BMD Longitud GPSurface
## 3.63439 0.26381 137.91720 37.09100 41.00512
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.69381 71.71914 99.70932 14.00000 2.00000Vectores de medias con 3 cluster.
colMeans(subset(Datos3_con_media[3:18], Cluster3 == 1))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 3.00000 0.00000 14.00000 306.80000 0.27303
## Tb.N BMC BMD Longitud GPSurface
## 3.63439 0.26381 137.91720 37.09100 41.00512
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.69381 71.71914 99.70932 14.00000 2.00000
## Cluster3
## 1.00000colMeans(subset(Datos3_con_media[3:18], Cluster3 == 2))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 3.0000000 230.7692308 14.0000000 216.6307692 0.2685615
## Tb.N BMC BMD Longitud GPSurface
## 3.6757538 0.1907077 122.8738462 34.1626923 34.4485692
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.5732538 71.7269413 99.1175006 14.0000000 1.0000000
## Cluster3
## 2.0000000colMeans(subset(Datos3_con_media[3:18], Cluster3 == 3))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 3.0000000 807.6923077 14.0000000 207.2153846 0.3732782
## Tb.N BMC BMD Longitud GPSurface
## 2.6334923 0.1679077 113.2255385 33.7961538 33.4378769
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.6814846 71.7625893 99.2139962 14.0000000 1.0000000
## Cluster3
## 3.0000000Vectores de medias con 4 cluster.
colMeans(subset(Datos3_con_media[3:19], Cluster4 == 1))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 3.000000 0.000000 14.000000 312.183333 0.251250
## Tb.N BMC BMD Longitud GPSurface
## 3.868433 0.268300 138.779333 37.171667 40.009517
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 4.855683 71.701187 99.725402 14.000000 2.000000
## Cluster3 Cluster4
## 1.000000 1.000000colMeans(subset(Datos3_con_media[3:19], Cluster4 == 2))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 3.0000000 230.7692308 14.0000000 216.6307692 0.2685615
## Tb.N BMC BMD Longitud GPSurface
## 3.6757538 0.1907077 122.8738462 34.1626923 34.4485692
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.5732538 71.7269413 99.1175006 14.0000000 1.0000000
## Cluster3 Cluster4
## 2.0000000 2.0000000colMeans(subset(Datos3_con_media[3:19], Cluster4 == 3))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 3.0000000 807.6923077 14.0000000 207.2153846 0.3732782
## Tb.N BMC BMD Longitud GPSurface
## 2.6334923 0.1679077 113.2255385 33.7961538 33.4378769
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.6814846 71.7625893 99.2139962 14.0000000 1.0000000
## Cluster3 Cluster4
## 3.0000000 3.0000000colMeans(subset(Datos3_con_media[3:19], Cluster4 == 4))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 3.000000 0.000000 14.000000 298.725000 0.305700
## Tb.N BMC BMD Longitud GPSurface
## 3.283325 0.257075 136.624000 36.970000 42.498525
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 6.951000 71.746067 99.685195 14.000000 2.000000
## Cluster3 Cluster4
## 1.000000 4.000000Conclusión número de clusters
Estudio 1
Decidimos elegir k=3. Por tanto, los datos se quedan agrupados en 3 clusters.
set.seed(123)
k3_1 <- kmeans(PCA_1scores_scale, centers = 3, nstart = 25)
fviz_cluster(k3_1, data = PCA_1scores_scale) + ggtitle("k = 3") 
Vector de medias de cada cluster:
colMeans(subset(Datos1_con_media[3:18], Cluster3 == 1))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.0000000 25.7600000 14.0000000 277.7183396 0.3170226
## Tb.N BMC BMD Longitud GPSurface
## 2.9973644 0.2509190 138.9594800 38.2376655 38.5942069
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.0396320 118.5800000 158.3110480 14.0000000 1.2000000
## Cluster3
## 1.0000000colMeans(subset(Datos1_con_media[3:18], Cluster3 == 2))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.0000000 36.4615385 14.0000000 303.2198839 0.2967179
## Tb.N BMC BMD Longitud GPSurface
## 3.3144923 0.2766495 145.1298718 39.0961538 41.5172458
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.7030430 122.4230769 171.1128739 14.0000000 1.7307692
## Cluster3
## 2.0000000colMeans(subset(Datos1_con_media[3:18], Cluster3 == 3))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 1.0000000 0.0000000 14.0000000 348.2800000 0.3444445
## Tb.N BMC BMD Longitud GPSurface
## 2.6883400 0.2810602 149.7147188 40.7620000 40.7213439
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 4.3462400 115.9400000 165.4812500 14.0000000 1.4000000
## Cluster3
## 3.0000000Estudio 2
Decidimos elegir k=4. Por tanto, los datos se quedan agrupados en 4 clusters.
set.seed(123)
k4_2 <- kmeans(PCA_2scores_scale, centers = 4, nstart = 25)
fviz_cluster(k4_2, data = PCA_2scores_scale) + ggtitle("k = 4") 
Vector de medias de cada cluster:
colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 1))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 8.8888889 21.0000000 341.0237037 0.3253313
## Tb.N BMC BMD Longitud GPSurface
## 2.7106444 0.3234977 176.6955556 41.7344344 45.6307988
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.9622333 120.6555556 161.0037037 28.0000000 2.0000000
## Cluster3 Cluster4
## 2.0000000 1.0000000colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 2))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 160.0000000 21.0000000 354.6500000 0.3189563
## Tb.N BMC BMD Longitud GPSurface
## 3.1786250 0.3242581 157.4050000 41.7466503 42.7150195
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.0384500 117.8583333 175.7666667 28.0000000 1.0000000
## Cluster3 Cluster4
## 1.0000000 2.0000000colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 3))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 70.0000000 21.0000000 321.5750000 0.3133542
## Tb.N BMC BMD Longitud GPSurface
## 2.8665125 0.3247514 151.8255000 41.7371055 41.3465789
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.1177875 111.9750000 151.6625000 28.0000000 1.8750000
## Cluster3 Cluster4
## 3.0000000 3.0000000colMeans(subset(Datos2.2_con_media[3:19], Cluster4 == 4))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 2.0000000 213.3333333 21.0000000 343.7444444 0.3089356
## Tb.N BMC BMD Longitud GPSurface
## 3.2059111 0.3250970 151.2567778 41.7493159 44.5133260
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 6.6399778 112.8333333 143.7666667 28.0000000 1.3333333
## Cluster3 Cluster4
## 2.1111111 4.0000000Estudio 3
Decidimos elegir k=3. Por tanto, los datos se quedan agrupados en 3 clusters.
set.seed(123)
k3_3 <- kmeans(PCA_3scores, centers = 3, nstart = 25)
fviz_cluster(k3_3, data = PCA_3scores) + ggtitle("k = 3")
Vector de medias de cada cluster:
colMeans(subset(Datos3_con_media[3:18], Cluster3 == 1))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 3.00000 0.00000 14.00000 306.80000 0.27303
## Tb.N BMC BMD Longitud GPSurface
## 3.63439 0.26381 137.91720 37.09100 41.00512
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.69381 71.71914 99.70932 14.00000 2.00000
## Cluster3
## 1.00000colMeans(subset(Datos3_con_media[3:18], Cluster3 == 2))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 3.0000000 230.7692308 14.0000000 216.6307692 0.2685615
## Tb.N BMC BMD Longitud GPSurface
## 3.6757538 0.1907077 122.8738462 34.1626923 34.4485692
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.5732538 71.7269413 99.1175006 14.0000000 1.0000000
## Cluster3
## 2.0000000colMeans(subset(Datos3_con_media[3:18], Cluster3 == 3))## Estudio Dosis DiasDosis PesoFinal Tb.Sp
## 3.0000000 807.6923077 14.0000000 207.2153846 0.3732782
## Tb.N BMC BMD Longitud GPSurface
## 2.6334923 0.1679077 113.2255385 33.7961538 33.4378769
## GPVolume PesoInicial PesoIntermedio DiasEspera Cluster2
## 5.6814846 71.7625893 99.2139962 14.0000000 1.0000000
## Cluster3
## 3.0000000ANÁLISIS DISCRIMINANTE LINEAL
ESTUDIO 1
set.seed(123)
PCA_1scores_scale$Cluster3 <- k3_1$cluster
datos<-PCA_1scores_scale[,-4]
datos$Cluster3 <- as.factor(datos$Cluster3)
set.seed(123)
trainIndex<-createDataPartition(datos$Cluster3,p=0.80)$Resample1
datos_train<-datos[trainIndex,]
datos_test<-datos[-trainIndex,]Nuestras variables predictoras van a ser las 3 primeras componentes principales y la dosis.
plot(datos_train[, 1:3], col = datos$Cluster3, pch = 19, main="Datos de entrenamiento")
plot(datos_test[, 1:3], col = datos$Cluster3, pch = 19, main="Datos de test")
En primer lugar exploramos cómo de bien (o mal) clasifica en los grupos de efectividad de la dosis cada una de las variables explicativas consideradas de forma independiente. Para ello dibujamos los histogramas superpuestos. Si los histogramas se separan, la variable considerada sería un buen clasificador individual para el grupo.
p1 <- ggplot(data = datos_train, aes(x = PC1, fill = Cluster3)) + geom_histogram(position = "identity", alpha = 0.5)
p2 <- ggplot(data = datos_train, aes(x = PC2, fill = Cluster3)) + geom_histogram(position = "identity", alpha = 0.5)
p3 <- ggplot(data = datos_train, aes(x = PC3, fill = Cluster3)) + geom_histogram(position = "identity", alpha = 0.5)
ggarrange(p1, p2, p3, nrow = 3, common.legend = TRUE)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Parece que la variable PC3 de la tercera componente principal es la que mejor separa entre los dos grupos.
Veamos qué pares de variables separan mejor:
pairs(x = datos_train[, c("PC1","PC2","PC3")],
col = c("blue", "orange")[datos_train$Cluster3], pch = 19)
Las combinaciones PC2-PC3 y PC1-PC3 parecen separar bien los datos.
Veamos como de bien separan las tres variables y si por tanto tiene sentido usarlas para hacer el análisis discriminante.
scatterplot3d(datos_train$PC1, datos_train$PC2, datos_train$PC3,
color = c("blue", "orange")[datos_train$Cluster3], pch = 19,
grid = TRUE, xlab = "var_exp1", ylab = "var_exp2",
zlab = "var_exp3", angle = 65, cex.axis = 0.6)
legend("topleft",
bty = "n", cex = 0.8,
title = "efecividad_dosis",
c("a", "b"), fill = c("blue", "orange"))
Separan bastante bien, así que tiene sentido proceder con el análisis discriminante.
A continuación hacemos una exploración gráfica de la normalidad de las distribuciones univariantes de nuestros predictores representando los histogramas y los gráficos qqplots.
par(mfcol = c(2, 3))
for (k in 1:3) {
j0 <- names(datos_train)[k]
x0 <- seq(min(datos_train[, k]), max(datos_train[, k]), le = 50)
for (i in 1:2) {
i0 <- levels(datos_train$Cluster3)[i]
x <- datos[datos_train$Cluster3 == i0, j0]
hist(x, proba = T, col = grey(0.8), main = paste("efectividad_dosis", i0), xlab = j0)
lines(x0, dnorm(x0, mean(x), sd(x)), col = "blue", lwd = 2)
}
}
##### Representation of normal quantiles of each variable for each species
par(mfrow=c(2,3))
for (k in 1:3) {
j0 <- names(datos_train)[k]
x0 <- seq(min(datos_train[, k]), max(datos_train[, k]), le = 50)
for (i in 1:2) {
i0 <- levels(datos_train$Cluster3)[i]
x <- datos[datos_train$Cluster3 == i0, j0]
qqnorm(x, main = paste("efectividad_dosis", i0, j0), pch = 19, col = i + 1)
qqline(x)
}
}
#####Normalidad univariante
datos_tidy <- melt(datos_train, value.name = "valor") ## Using Cluster3 as id variablesdatos_tidy %>%
group_by(Cluster3,variable) %>%
summarise(p_value_Shapiro.test = round(shapiro.test(valor)$p.value,5))## `summarise()` has grouped output by 'Cluster3'. You can override using the
## `.groups` argument.## # A tibble: 9 × 3
## # Groups: Cluster3 [3]
## Cluster3 variable p_value_Shapiro.test
## <fct> <fct> <dbl>
## 1 1 PC1 0.590
## 2 1 PC2 0.0597
## 3 1 PC3 0.136
## 4 2 PC1 0.700
## 5 2 PC2 0.370
## 6 2 PC3 0.438
## 7 3 PC1 0.243
## 8 3 PC2 0.524
## 9 3 PC3 0.284Observamos que el p-valor>0.05, podemos asumir normalidad univariante en cada variable.
Tests de normalidad multivariante:
outliers <- mvn(data = datos_train[,-4], mvnTest = "hz", multivariateOutlierMethod = "quan")
royston_test <- mvn(data = datos_train[,-4], mvnTest = "royston", multivariatePlot = "qq")
royston_test$multivariateNormality## Test H p value MVN
## 1 Royston 8.019127 0.04562044 NOhz_test <- mvn(data = datos_train[,-4], mvnTest = "hz")
hz_test$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 0.6837849 0.3736335 YESHomogeneidad de varianzas
boxM(data = datos_train[,-4], grouping = datos_train[, 4])##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: datos_train[, -4]
## Chi-Sq (approx.) = 13.061, df = 12, p-value = 0.3646Observamos que se da normalidad y homogeneidad de las varianzas.
modelo_lda <- lda(formula = Cluster3 ~ PC1 + PC2 + PC3, data=datos_train)
modelo_lda## Call:
## lda(Cluster3 ~ PC1 + PC2 + PC3, data = datos_train)
##
## Prior probabilities of groups:
## 1 2 3
## 0.08888889 0.46666667 0.44444444
##
## Group means:
## PC1 PC2 PC3
## 1 -1.3718383 -0.4563999 -0.3148362
## 2 -0.5783738 0.3208546 0.4149847
## 3 0.8157435 -0.1246100 -0.4954859
##
## Coefficients of linear discriminants:
## LD1 LD2
## PC1 1.6060635 -0.4559855
## PC2 -0.6357729 -0.6652696
## PC3 -0.7826050 -0.6374750
##
## Proportion of trace:
## LD1 LD2
## 0.9513 0.0487#probabilidades: 0.6198, 0.3902Una vez construido el clasificador, podemos clasificar nuevos datos en función de sus medidas sin más que llamar a la función predict.
Veamos para todos los datos test.
nuevas_observaciones<-datos_test
predict(object = modelo_lda, newdata = nuevas_observaciones) ## $class
## [1] 2 3 2 2 2 2 2 3 3 3 3
## Levels: 1 2 3
##
## $posterior
## 1 2 3
## 4 6.372456e-02 9.347727e-01 0.0015027191
## 7 1.640030e-02 1.226638e-01 0.8609359134
## 13 1.605429e-02 9.835519e-01 0.0003938489
## 21 1.960784e-02 9.606043e-01 0.0197878602
## 24 2.714848e-02 9.721102e-01 0.0007413120
## 25 2.748075e-01 7.146342e-01 0.0105582776
## 33 4.689169e-02 7.266001e-01 0.2265082433
## 37 1.281882e-05 8.759697e-05 0.9998995842
## 42 3.187139e-04 4.470158e-03 0.9952111282
## 50 1.124855e-05 7.747008e-04 0.9992140507
## 52 6.855588e-04 1.749723e-03 0.9975647187
##
## $x
## LD1 LD2
## 4 -1.7567557 -0.3552411
## 7 0.7677199 0.5384874
## 13 -2.1031712 -1.4725951
## 21 -0.9112295 -1.1209042
## 24 -1.9357527 -1.0470107
## 25 -1.1737072 1.0225313
## 33 -0.1440404 -0.1452081
## 37 3.0223824 0.9563182
## 42 1.8586148 0.2403761
## 50 2.4826781 -0.8464948
## 52 2.0525918 1.5382039Validación del modelo
pred <- predict(object = modelo_lda, newdata = datos_test)
confusionmatrix(datos_test$Cluster3, pred$class)## new 1 new 2 new 3
## 1 0 1 0
## 2 0 5 0
## 3 0 0 5Clasificación del porcentaje del error
trainig_error <- mean(datos_test$Cluster3 != pred$class) * 100
paste("trainig_error=", trainig_error, "%")## [1] "trainig_error= 9.09090909090909 %"Visualización de las clasificaciones
partimat(Cluster3 ~ PC1 + PC2 + PC3,
data = datos_test, method = "lda", prec = 200,
image.colors = c("green", "orange","blue"),
col.mean = "yellow",nplots.vert =1, nplots.hor=3)
Training error: 9%. Cierto error en PC1-PC2, 18%; PC1-PC3, 9%.
ESTUDIO 2
PCA_2scores_scale$Cluster4 <- k4_2$cluster
datos<-PCA_2scores_scale[,-4]
datos$Cluster4 <- as.factor(datos$Cluster4)
# Partitioning the dataset: training (80%) + test (20%)
set.seed(123)
trainIndex<-createDataPartition(datos$Cluster4,p=0.80)$Resample1
datos_train<-datos[trainIndex,]
datos_test<-datos[-trainIndex,]Nuestras variables predictoras van a ser las 3 primeras componentes principales y la dosis.
plot(datos_train[, 1:3], col = datos$Cluster3, pch = 19, main="Datos de entrenamiento")
plot(datos_test[, 1:3], col = datos$Cluster3, pch = 19, main="Datos de test")
Exploramos clasificación.
p1 <- ggplot(data = datos_train, aes(x = PC1, fill = Cluster4)) + geom_histogram(position = "identity", alpha = 0.5)
p2 <- ggplot(data = datos_train, aes(x = PC2, fill = Cluster4)) + geom_histogram(position = "identity", alpha = 0.5)
p3 <- ggplot(data = datos_train, aes(x = PC3, fill = Cluster4)) + geom_histogram(position = "identity", alpha = 0.5)
ggarrange(p1, p2, p3, nrow = 3, common.legend = TRUE)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Parece que la variable PC3 de la tercera componente principal es la que mejor separa entre los dos grupos.
Veamos qué pares de variables separan mejor:
pairs(x = datos_train[, c("PC1","PC2","PC3")],
col = c("blue", "orange")[datos_train$Cluster4], pch = 19)
Las combinaciones PC2-PC3 y PC1-PC3 parecen separar bien los datos
Veamos como de bien separan las tres variables y si por tanto tiene sentido usarlas para hacer el análisis discriminante
scatterplot3d(datos_train$PC1, datos_train$PC2, datos_train$PC3,
color = c("blue", "orange")[datos_train$Cluster4], pch = 19,
grid = TRUE, xlab = "var_exp1", ylab = "var_exp2",
zlab = "var_exp3", angle = 65, cex.axis = 0.6)
legend("topleft",
bty = "n", cex = 0.8,
title = "efecividad_dosis",
c("a", "b"), fill = c("blue", "orange"))
Separan bastante bien, así que tiene sentido proceder con el análisis discriminante.
A continuación hacemos una exploración gráfica de la normalidad de las distribuciones univariantes de nuestros predictores representando los histogramas y los gráficos qqplots.
par(mfcol = c(2, 3))
for (k in 1:3) {
j0 <- names(datos_train)[k]
x0 <- seq(min(datos_train[, k]), max(datos_train[, k]), le = 50)
for (i in 1:2) {
i0 <- levels(datos_train$Cluster4)[i]
x <- datos[datos_train$Cluster4 == i0, j0]
hist(x, proba = T, col = grey(0.8), main = paste("efectividad_dosis", i0), xlab = j0)
lines(x0, dnorm(x0, mean(x), sd(x)), col = "blue", lwd = 2)
}
}
par(mfrow=c(2,3))
for (k in 1:3) {
j0 <- names(datos_train)[k]
x0 <- seq(min(datos_train[, k]), max(datos_train[, k]), le = 50)
for (i in 1:2) {
i0 <- levels(datos_train$Cluster4)[i]
x <- datos[datos_train$Cluster4 == i0, j0]
qqnorm(x, main = paste("efectividad_dosis", i0, j0), pch = 19, col = i + 1)
qqline(x)
}
}
#Normalidad univariante
datos_tidy <- melt(datos_train, value.name = "valor") ## Using Cluster4 as id variablesdatos_tidy %>%
group_by(Cluster4,variable) %>%
summarise(p_value_Shapiro.test = round(shapiro.test(valor)$p.value,5))## `summarise()` has grouped output by 'Cluster4'. You can override using the
## `.groups` argument.## # A tibble: 12 × 3
## # Groups: Cluster4 [4]
## Cluster4 variable p_value_Shapiro.test
## <fct> <fct> <dbl>
## 1 1 PC1 0.349
## 2 1 PC2 0.641
## 3 1 PC3 0.586
## 4 2 PC1 0.146
## 5 2 PC2 0.149
## 6 2 PC3 0.766
## 7 3 PC1 0.530
## 8 3 PC2 0.562
## 9 3 PC3 0.795
## 10 4 PC1 0.813
## 11 4 PC2 0.792
## 12 4 PC3 0.367Observamos que el p-valor>0.05, podemos asumir normalidad univariante en cada.
Tests de normalidad multivariante:
outliers <- mvn(data = datos_train[,-4], mvnTest = "hz", multivariateOutlierMethod = "quan")
royston_test <- mvn(data = datos_train[,-4], mvnTest = "royston", multivariatePlot = "qq")
royston_test$multivariateNormality## Test H p value MVN
## 1 Royston 1.998275 0.5727658 YEShz_test <- mvn(data = datos_train[,-4], mvnTest = "hz")
hz_test$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 0.6472712 0.3817672 YESHomogeneidad de varianzas
boxM(data = datos_train[,-4], grouping = datos_train[, 4])##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: datos_train[, -4]
## Chi-Sq (approx.) = 20.802, df = 18, p-value = 0.2895Observamos que se da normalidad y homogeneidad de las varianzas.
modelo_lda <- lda(formula = Cluster4 ~ PC1 + PC2 + PC3, data=datos_train)
modelo_lda## Call:
## lda(Cluster4 ~ PC1 + PC2 + PC3, data = datos_train)
##
## Prior probabilities of groups:
## 1 2 3 4
## 0.21875 0.25000 0.25000 0.28125
##
## Group means:
## PC1 PC2 PC3
## 1 0.6599792 0.03286496 1.0466227
## 2 -0.5456697 0.91151269 0.4214908
## 3 -1.1064918 -1.02213585 -0.2658523
## 4 0.5758493 0.03657006 -1.0758341
##
## Coefficients of linear discriminants:
## LD1 LD2 LD3
## PC1 -1.3886224 0.97816263 -0.5633177
## PC2 -1.0444490 -0.04111225 0.9622804
## PC3 -0.8087304 -1.34473390 -0.4414437
##
## Proportion of trace:
## LD1 LD2 LD3
## 0.5642 0.3511 0.0847#probabilidades: 0.9513, 0.0487Una vez construido el clasificador, podemos clasificar nuevos datos en función de sus medidas sin más que llamar a la función predict
Veamos para todos los datos test
nuevas_observaciones<-datos_test
predict(object = modelo_lda, newdata = nuevas_observaciones) ## $class
## [1] 1 4 4 3 2 2
## Levels: 1 2 3 4
##
## $posterior
## 1 2 3 4
## 2 0.9827994437 0.0171845072 5.716587e-06 1.033252e-05
## 11 0.3769804388 0.0025946913 8.099923e-05 6.203439e-01
## 19 0.0055702677 0.0003012557 2.306756e-07 9.941282e-01
## 21 0.0003815275 0.1069681004 8.917523e-01 8.980249e-04
## 33 0.3384897107 0.4389712756 7.762480e-07 2.225382e-01
## 38 0.2714866078 0.4313807576 4.669169e-03 2.924635e-01
##
## $x
## LD1 LD2 LD3
## 2 -2.1123001 -2.2635332 -2.04196261
## 11 -1.3713408 1.6426389 -2.26599200
## 19 -1.9791146 3.0688773 -0.39825960
## 21 1.8261149 -1.4142720 0.75168656
## 33 -2.7972879 0.9137039 1.49106874
## 38 -0.6681194 0.2290076 -0.03096111Validación del modelo
pred <- predict(object = modelo_lda, newdata = datos_test)
confusionmatrix(datos_test$Cluster4, pred$class)## new 1 new 2 new 3 new 4
## 1 1 0 0 0
## 2 0 2 0 0
## 3 0 0 1 0
## 4 0 0 0 2Clasificación del porcentaje del error
trainig_error <- mean(datos_test$Cluster4 != pred$class) * 100
paste("trainig_error=", trainig_error, "%")## [1] "trainig_error= 0 %"Visualización de las clasificaciones
partimat(Cluster4 ~ PC1 + PC2 + PC3,
data = datos_test, method = "lda", prec = 200,
image.colors = c("green", "orange","blue","red"),
col.mean = "yellow",nplots.vert =1, nplots.hor=3)
Training error: 0%. Rrror de 40% en PC1-PC3, sustancial por falta de datos.
ESTUDIO 3
PCA_3scores$Cluster3 <- k3_3$cluster
datos<-PCA_3scores
datos$Cluster3 <- as.factor(datos$Cluster3)
# Partitioning the dataset: training (80%) + test (20%)
set.seed(123)
trainIndex<-createDataPartition(datos$Cluster3,p=0.80)$Resample1
datos_train<-datos[trainIndex,]
datos_test<-datos[-trainIndex,]Nuestras variables predictoras van a ser las 3 primeras componentes principales y la dosis.
plot(datos_train[, 1:3], col = datos$Cluster3, pch = 19, main="Datos de entrenamiento")
plot(datos_test[, 1:3], col = datos$Cluster3, pch = 19, main="Datos de test")
Exploramos clasificación.
p1 <- ggplot(data = datos_train, aes(x = PC1, fill = Cluster3)) + geom_histogram(position = "identity", alpha = 0.5)
p2 <- ggplot(data = datos_train, aes(x = PC2, fill = Cluster3)) + geom_histogram(position = "identity", alpha = 0.5)
p3 <- ggplot(data = datos_train, aes(x = PC3, fill = Cluster3)) + geom_histogram(position = "identity", alpha = 0.5)
ggarrange(p1, p2, p3, nrow = 3, common.legend = TRUE)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Parece que la variable PC3 de la tercera componente principal es la que mejor separa entre los dos grupos.
Veamos qué pares de variables separan mejor:
pairs(x = datos_train[, c("PC1","PC2","PC3")],
col = c("blue", "orange")[datos_train$Cluster3], pch = 19)
Las combinaciones PC2-PC3 y PC1-PC3 parecen separar bien los datos
Veamos como de bien separan las tres variables y si por tanto tiene sentido usarlas para hacer el análisis discriminante
scatterplot3d(datos_train$PC1, datos_train$PC2, datos_train$PC3,
color = c("blue", "orange")[datos_train$Cluster3], pch = 19,
grid = TRUE, xlab = "var_exp1", ylab = "var_exp2",
zlab = "var_exp3", angle = 65, cex.axis = 0.6)
legend("topleft",
bty = "n", cex = 0.8,
title = "efecividad_dosis",
c("a", "b"), fill = c("blue", "orange"))
Separan bastante bien, así que tiene sentido proceder con el análisis discriminante.
A continuación hacemos una exploración gráfica de la normalidad de las distribuciones univariantes de nuestros predictores representando los histogramas y los gráficos qqplots.
par(mfcol = c(2, 3))
for (k in 1:3) {
j0 <- names(datos_train)[k]
x0 <- seq(min(datos_train[, k]), max(datos_train[, k]), le = 50)
for (i in 1:2) {
i0 <- levels(datos_train$Cluster3)[i]
x <- datos[datos_train$Cluster3 == i0, j0]
hist(x, proba = T, col = grey(0.8), main = paste("efectividad_dosis", i0), xlab = j0)
lines(x0, dnorm(x0, mean(x), sd(x)), col = "blue", lwd = 2)
}
}
par(mfrow=c(2,3))
for (k in 1:3) {
j0 <- names(datos_train)[k]
x0 <- seq(min(datos_train[, k]), max(datos_train[, k]), le = 50)
for (i in 1:2) {
i0 <- levels(datos_train$Cluster3)[i]
x <- datos[datos_train$Cluster3 == i0, j0]
qqnorm(x, main = paste("efectividad_dosis", i0, j0), pch = 19, col = i + 1)
qqline(x)
}
}
#Normalidad univariante
datos_tidy <- melt(datos_train, value.name = "valor") ## Using Cluster3 as id variablesdatos_tidy %>%
group_by(Cluster3,variable) %>%
summarise(p_value_Shapiro.test = round(shapiro.test(valor)$p.value,5))## `summarise()` has grouped output by 'Cluster3'. You can override using the
## `.groups` argument.## # A tibble: 9 × 3
## # Groups: Cluster3 [3]
## Cluster3 variable p_value_Shapiro.test
## <fct> <fct> <dbl>
## 1 1 PC1 0.110
## 2 1 PC2 0.348
## 3 1 PC3 0.0920
## 4 2 PC1 0.0484
## 5 2 PC2 0.0903
## 6 2 PC3 0.257
## 7 3 PC1 0.217
## 8 3 PC2 0.799
## 9 3 PC3 0.212Observamos que el p-valor>0.05, podemos asumir normalidad univariante en todas excepto en una.
Tests de normalidad multivariante:
outliers <- mvn(data = datos_train[,-4], mvnTest = "hz", multivariateOutlierMethod = "quan")
royston_test <- mvn(data = datos_train[,-4], mvnTest = "royston", multivariatePlot = "qq")
royston_test$multivariateNormality## Test H p value MVN
## 1 Royston 9.794117 0.02039988 NOhz_test <- mvn(data = datos_train[,-4], mvnTest = "hz")
hz_test$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 1.175856 0.001754329 NOHomogeneidad de varianzas
boxM(data = datos_train[,-4], grouping = datos_train[, 4])##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: datos_train[, -4]
## Chi-Sq (approx.) = 16.557, df = 12, p-value = 0.167Observamos que se da normalidad y homogeneidad de las varianzas.
modelo_lda <- lda(formula = Cluster3 ~ PC1 + PC2 + PC3, data=datos_train)
modelo_lda## Call:
## lda(Cluster3 ~ PC1 + PC2 + PC3, data = datos_train)
##
## Prior probabilities of groups:
## 1 2 3
## 0.3666667 0.2666667 0.3666667
##
## Group means:
## PC1 PC2 PC3
## 1 2.0553136 -0.6432026 -0.3127310
## 2 -3.4443121 -0.4688611 -0.1386789
## 3 0.3644586 1.0726328 0.6388531
##
## Coefficients of linear discriminants:
## LD1 LD2
## PC1 1.2431245 0.05057212
## PC2 0.1505819 -1.07116237
## PC3 0.1192883 -0.72051310
##
## Proportion of trace:
## LD1 LD2
## 0.8462 0.1538#probabilidades:Una vez construido el clasificador, podemos clasificar nuevos datos en función de sus medidas sin más que llamar a la función predict
Veamos para todos los datos test
nuevas_observaciones<-datos_test
predict(object = modelo_lda, newdata = nuevas_observaciones) ## $class
## [1] 3 3 2 2 1 1
## Levels: 1 2 3
##
## $posterior
## 1 2 3
## 4 2.139487e-02 1.761181e-07 9.786049e-01
## 7 1.632248e-02 7.016742e-11 9.836775e-01
## 10 4.180382e-10 9.999998e-01 1.807804e-07
## 16 2.524699e-10 9.999995e-01 4.887665e-07
## 30 9.999656e-01 1.401340e-15 3.441981e-05
## 32 9.997673e-01 3.234498e-11 2.326927e-04
##
## $x
## LD1 LD2
## 4 0.9036463 -1.2073013
## 7 2.0945594 -2.1026849
## 10 -4.2218932 1.3303718
## 16 -4.2477080 0.7730017
## 30 3.9206753 2.1962734
## 32 2.4179979 2.4601340Validación del modelo
pred <- predict(object = modelo_lda, newdata = datos_test)
confusionmatrix(datos_test$Cluster3, pred$class)## new 1 new 2 new 3
## 1 2 0 0
## 2 0 2 0
## 3 0 0 2Clasificación del porcentaje del error
trainig_error <- mean(datos_test$Cluster3 != pred$class) * 100
paste("trainig_error=", trainig_error, "%")## [1] "trainig_error= 0 %"Visualización de las clasificaciones
partimat(Cluster3 ~ PC1 + PC2 + PC3,
data = datos_test, method = "lda", prec = 200,
image.colors = c("green", "orange","blue"),
col.mean = "yellow",nplots.vert =1, nplots.hor=3)
Training error: 0%. Error en PC2-PC3 del 33’3%.